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Introduction 

Description 

The Numerical Analysis Library provides you with quick access to 57 commonly used numer- 
ical analysis routines. These routines are in subprogram form. They can be called up from the 
library as you need them and appended to your program in memory. Except for error mes- 
sages, these subprograms contain no input or output operations. Drivers are provided with 
each subprogram to illustrate some of the techniques for calling the subprogram. 

Each section of this manual contains several techniques, giving you the flexibility to choose 
the most appropriate routine for your specific need. 

System Configuration 

To use the Numerical Analysis Library, you need: 

• a 9826 or 9836 computer with a BASIC language system. 

• Numerical Analysis software — part number 98821A. This includes a user's manual and 
two flexible discs containing the programs. 

Using the Numerical Analysis Routines 

Drivers and subprograms are contained in mass storage files with the same names. Driver file 
names have all capital letters and subprogram file names are initially capped followed by 
lower case letters. For example, "MULLER" contains the driver for subprogram Muller con- 
tained in file "Muller". 

To get a particular numerical analysis routine, insert the appropriate Numerical Analysis flexi- 
ble disc into the disc drive with the machine turned on and the BASIC language system 
loaded. Either select a suitable routine, turn to its user instructions in this manual, and follow 
them, or write your own "driver" program, including the statement LOADSUB ALL FROM 
"<filename>" where <filename> is the appropriate numerical analysis subprogram to be 
used by your program and execute it. If you accidently "loadsub" the wrong file and if that 
file has any subprograms (either SUB or DEF FN segments), those subprograms will be 
loaded, right or wrong. If you "loadsub" a file that has been "SAVE'd" (or "RE-SAVE' d") as 
opposed to "STORE'd" (or "RE-STORE' d") you will get an error 58. If you "loadsub" a file 
that contains no subprograms, no error will appear and nothing will be loaded. 

Using the Numerical Analysis Drivers 

The driver programs allow you to use the routines without writing a main program. To use 
those drivers, type: LOAD "<file name containing driver>", press EXECUTE and RUN. 
Necessary subprograms will be automatically "loaded", and the program will ask you for your 
entries. 
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Redefining User Functions 

Some programs (Root Finders, Integration, Differential Equations) need a user function (con- 
tained in files "KRYWEN" or "FUNC".) You have to STORE those files before "loading" 
and "running" the driver. However, if you want to check or redefine your function after you 
have LOADed the driver and started to RUN it: 

1. PAUSE the current program. 

2. Type EDIT F and press EXECUTE. 

3. EDIT mode displays the program line containing your function. 

4. Press ENTER after modifying the program line. 

5. Press RUN to restart the program. 
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Introduction 

Description 

This section contains four different algorithms for finding roots of equations. Unfortunately, 
solving nonlinear equations is very difficult in practice. The methods given here are local 
methods. In the presence of good starting guesses, these work quite well. However, in the 
absence of good estimates, these techniques range from wild divergence to rapid convergence 
depending upon the particular user-defined function. 

Hence, a word of caution. There is no such thing as a fool-proof nonlinear technique. When 
possible, graphic or analytic means should be used to get a rough estimate of the roots. In 
terms of speed, stability and precision, each of the above methods has its own advantages 
and disadvantages. With a number of different techniques to choose from, you may find one 
that will most adequately meet your needs. 

Routines 

• Bisect - an iterative rootfinder using the bisection method. 

• Roota - a general rootfinder using a modified secant method. 

• Work - finds the root of a user-defined real-valued continuous function by a modified 
secant method; contained in file "Roota". 

• Chkit - checks if a root has been found or bracketed; contained in file "Roota". 

• Monton - checks the three points, (x x , yj, (x 2 , y 2 ) and (x 3 , y 3 ), to see that y = f(x) is 
strictly increasing or decreasing; contained in file "Roota". 

• FN Incres - used to check if too many marching steps have been taken in searching for a 
change of sign; contained in file "Roota". 

• Muller - a general rootfinder using a quadratic method. 

• Siljak - a polynomial rootfinder with complex coefficients. 
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Bisection Method 

Description 

This subprogram will search for solutions of f(x) = over an interval [a,b] where you define 
the continuous real-valued function f(x). The function may be algebraic of the form 



(a + ajx 8 ' + a 2 x e2 + +a n x e ") with a; real and e, rational (e.g., x 

scendental (e.g., sin(x) + cos (x)). 



+ 3x' ' + x) or tran- 



You must specify the initial domain value, search increment, error tolerance, and maximum 
interval-halvings to be used. 

Program Usage 

Driver Utilization 

• File Name - "BISECT", disc 1 

The driver "BISECT" sets up the necessary input parameters for the subprogram Bisect and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Bisect", disc 1 

• Calling Syntax 

CALL Biset (A, B, Maxbi, Tol, Deltax,Root(*), F_(*), Err(*), Nroots) 

• Input Parameters 



A 

B 

Maxbi 

Tol 

Deltax 

Nroots 
• Output Parameters 
Root(*) 

F_(*) 
Err(*) 



Lower bound of search interval. 

Upper bound of search interval. 

Maximum number of bisections for each subinterval. 

Error tolerance for a root; a root x is accepted if 
|f(x)| < Tol*(l MAX x). 

Search increment; the subprogram begins at the lower 
bound and compares functional values at the ends of each 
subinterval, [a + iAx, a + (i + l)Ax], for a change of sign. 
Number of roots to be found. 



Vector containing the roots of the function; the value 
9.999999E99 is supplied for any roots not located. 

Vector containing the functional evaluations of the roots; 
again, the value 9.999999E99 is supplied for any roots not 
found. 

Vector containing the estimated accuracy of the roots; the 
value 9.999999E99 is supplied for any roots not found. 
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• Subprograms Required 

FN F(X) Given a domain value x, this function subprogram should 

return the range value of the function you supply. 



Special Considerations and Programming Hints 

• Because the number of roots located cannot be computed ahead of time, all roots are 
initialized to 9.999999E99 at the beginning of the subprogram. Thus, on output, any 
roots not found will contain the value 9.999999E99. 

• One of the advantages of this method is that whenever a root is found, both the accuracy 
of the root (contained in array Err(*)) and the accuracy of the functional evaluation of 
the root (contained in array F(*)) are known precisely. 

• Clearly, with enough effort, one can always locate a root to any desired accuracy with 
this subprogram. But compared with some of the other methods contained in this sec- 
tion, the bisection method converges rather slowly since it does not make full use of the 
information about f(x) available at each step. So, when possible, this method should not 
be used if there are a large number of roots to be found or if the rootfinder subprogram 
is to be called a number of times. 

• To define the function f(x): 

1. PAUSE the current program (if one is running). 

2. Type SCRATCH and press EXECUTE. 

3. In EDIT mode, type in the following function subprogram pressing ENTER after 
each line: 

10 DEF FNF(X) 

20 F: F = <user-defined function of the form ( X- 3)*(X + 4)> 

30 RETURN F 

40 FNEND 

4. RE-STORE the function subprogram on file "KRYWEN". 

5. LOAD file "BISECT" and press RUN. 

(See "Redefining User Functions" in the Introduction, following the Table of Contents, for 
another way to modify the function f(x).) 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Bisect." 
A = B = Maxbi = 

Tol = Deltax = Nroots = 

You may correct the data from the keyboard (e.g., Tol = 1E-6, EXECUTE). When 
CONTINUE is pressed, the program will resume execution at the next line. 
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• If, in a particular subinterval, the maximum number of iterations is exceeded, the sub- 
program will print the following warning message and the approximate root and accu- 
racies so far obtained: 

"MAX # OF BISECTIONS REACHED ON ROOT #<n>". 
"X BETWEEN <left endpoint> AND <right endpoint>" 
"F(X) = <y>" 

"ACCURACY TO (<rightendpoint> - <left endpoint>)" 
"AVERAGE VALUE STORED IN Root(*) AS APPROXIMATE X" 

The subprogram will then move to the next search interval. 

Methods and Formulae 

The subprogram Bisect will search for solutions of f(x) — over an interval [a,b] where you 
define the real-valued continuous function f(x). The function may be algebraic of the form 
a + ajx 6 ' + a 2 x 62 H + a n x Sn with a t real and e x rational (e.g., x J + 3x 3/2 + x) or tran- 
scendental (e.g., sin(x) + cos(x)). 

You must specify the search increment Ax and the error tolerance for f(x). The program then 
begins at the left of the interval and compares functional values at the ends of the subinterval 
[a, a + Ax]. If the functional values are of opposite sign then the bisection method is used to 
locate the root. Each subinterval [a + iAx, a + (i + l)Ax] is examined for a possible root. At 
most, one root per interval will be located and if there are multiple roots per interval, none 
will be located. You must also specify a maximum number of interval-halvings or bisections 
(Maxbi), so that an error tolerance that is not satisfied will result in the root localized to an 
interval of size 2~ Maxbl (b -a). The subprogram will examine N = INT b - a intervals. 

Ax 



Reference 

1. Stark, Peter A., Introduction to Numerical Methods, London: MacMillan Company, 
Collier-MacMillan Limited, 1970, pp. 95-96. 
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User Instructions 

1. If you have not already defined the function you wish to solve, you should do so at this 
time. 

2. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "BISECT" and press RUN. 

3. You will be asked to supply entries for the following items: 

- number of roots to be found 

- lower bound of search interval 

- upper bound of search interval 

- maximum number of bisections for searching for any one root in a subinterval 

- error tolerance where |f(x)|<Tol*(l MAX x) is the tolerance criterion 

- search increment or Ax. 

Press CONTINUE after each entry. 

4. The roots x, their functional values f(x), and the accuracy to which the root is found will 
be printed for all roots found over the interval. The value 9.999999E99 is inserted for 
any roots not found. 

Example 

User-defined function: 



User entries: 



Results: 



SIN(X) - COS(X)/(l + X*X) 



••!!■• OF ROOTS TO BE FOUND™ 4 

LOWER BOUND™ 

UPPER BOUND™ 10 

MAXIMUM * OF BISECTIONS** 20 

ERROR TOLERANCE™ i . E-8 

SEARCH INCREMENT™ . i 



ROOT FUNCTION VALUE ACCURACY 

6.23899&E-0i -9 . 442352E ■■••()? 1 . 9 7349E-0 7 

3 . 228892E+Q0 i . 620545E-08 1 . 907349E-07 

6 . 307698E+00 ■•■•! . 632362E--03 .1. . 907349E ™07 

9 . 43S884E+0 8 . 9S8369E-Q8 3. 051758E-06 
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Modified Secant Method 

Description 

Given a first guess, this subprogram will search for a solution of f(x) : 
continuous real-valued function f(x). 



where you define the 



Roots are found by marching along at a given step size until a change of sign is encountered. 
Then a modified secant method is employed to determine the zero of the function. 

You are required to specify the error tolerances for the root and for the function evaluation, as 
well as the step size and the maximum number of steps and iterations allowed. 

Program Usage 

Driver Utilization 

• File Name - "ROOTA", disc 1 

The driver "ROOTA" sets up the necessary input parameters for the subprogram Roota and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Roota", disc 1 

• Calling Syntax 

CALL Roota (Root, Err, Frsges, Step, Istpmx, Itrmx, Tola, Tolf) 

• Input Parameters 



Frsges 
Step 

Istpmx 

Itrmx 

Tola 
Tolf 



Output Parameters 

Root 
Err 



Initial guess for the root. 

Step size for marching when looking for an interval in which 
the function changes sign; step size must be positive. 

Maximum number of allowed steps in marching to find an 
interval in which the function changes sign. 

Maximum number of iterations allowed in subprogram work 
after an interval has been found in which the function 
changes sign. 

Tolerance for the root; i.e., if x : and x 2 are two consecutive 
approximations for the root, the average of Xj and x 2 is 
accepted as a root if |x t - x 2 | *£ (1 MAX (|xj| MAX |x 2 |)) * 
Tola. 

Tolerance for the function; i.e., an approximation x is 
accepted as a root if |f(x)| «£ Tolf. 



Root of f(x) 
f(Root) 
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• Subprograms Required 

Chkit Checks if a root is found or bracketed. 

FN F(X) User-defined function which, when given a domain value x, 

returns the value f(x). 

FN Incres Checks if too many marching steps have been taken in sub- 

program Roota. 

Monton Checks that the function is well-behaved in the search in- 

terval. 

Work Finds the root of the user-defined real-valued continuous 

function f by a modified secant method; requires that the 
functional values of the endpoints of the interval be of 
opposite sign. 

Further explanation of these subprograms is given on pages 1-8 and 1-10 through 1-15. 

Special Considerations and Programming Hints 

• Upon entry into subprogram Roota, there is a bad data check. If the check detects any 
"nonsense" data, an error message will be printed and the program pauses: 

"ERROR IN SUBPROGRAM Roota." 
Step = Istpmx = Itrmx = 

Tola = Tolf = 

The data may be corrected from the keyboard (e.g., Tola = IE -6, EXECUTE). When 
CONTINUE is pressed, the program will resume execution at the next line. 

• This method is slower than subprogram Muller, but has the advantage of having less 
trouble with functions that have somewhat unusual behavior, or functions that have 
periodic oscillations. 

• Subprogram Monton is used to check that the function is monotonic in the interval under 
consideration. If a change in direction is encountered, the following error message is 
printed and the program pauses: 

"ERROR IN SUBPROGRAM Monton." 

"Y = F(X) IS NOT MONOTONIC OR IS VERY FLAT." 

"SUGGEST TRYING SMALLER STEP SIZE OR DIFFERENT FIRST GUESS." 

XI = X2 = X3 = 

Yl = Y2 = Y3 = 

There are many possible causes of this error, the most common being that (X2, Y2) is a 
relative minimum or maximum. Choosing XI or X3 as the first guess and a step size less 
than |X1 - X2| may solve this problem. You may want to restart the program and make 
these changes. 

If the data indicates that the function is very flat in the interval, an initial guess on either 
side of the interval may help. 
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• To define the function f(x): 

1. PAUSE the current program (if one is running). 

2. Type SCRATCH and press EXECUTE. 

3. In EDIT mode, type in the following function subprogram pressing ENTER after 
each line: 

10 DEF FNF(X) 

20 F: F = < user-defined function of the form (X-3)*(X + 4) > 

30 RETURN F 

40 FNEND 

4. RE-STORE the function subprogram on file "KRYWEN". 

5. LOAD file "ROOTA" and press RUN. 

(See "Redefining User Functions" in the Introduction, following the Table of Contents, for 
another way to modify the function f(x). ) 

Methods and Formulae 

Subprogram Roota begins by initializing X2 - Frsges. Then subprogram Chkit is called. In 
Chkit a new x is generated and |f(x)| =£ Tolf is tested. If it is true, then x is accepted as the 
root. Otherwise a step to the left and a step to the right are tested. If a root is not found and if 
f (x) does not change sign, then a test is made to see which way to march: 

(y 2 - yj* yi < 0, then march right 
(y 2 - yi)* yi> 0, then march left. 

If a root or a change of sign is not found after Istpmx steps, an error message is printed and 
the program pauses. 

If a change of sign is found, then Chkit calls subprogram Work, where Root is found by a 
modified secant method which requires 2 ordinates to be of opposite sign. 

x 2 = (y 3 * xj - yi * x 3 )/(y 3 - y a ) 

A root is accepted if 

|f(x 2 )| s£ Tolf 
or 

|x x - x 2 | *£ (1 MAX (|x x | MAX |x 2 |)) * Tola 
where x 1 , x 2 are consecutive attempts at finding a root. 

Subprogram Monton is invoked repeatedly to test that f(x) is strictly increasing or decreasing 
and that x^ x 2 , x 3 are in ascending order. 
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User Instructions 

1. If you have not already defined the function you wish to solve, you should do so at this 
time. 

2. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "ROOTA" and press RUN. 

3. You will be asked to supply entries for the following items: 

- first guess for the root 

- step size (a positive value) for marching when looking for an interval in which the 
function changes sign 

- maximum number of steps allowed in marching to find an interval in which the 
function changes sign 

- maximum number of iterations allowed after an interval has been found in which the 
function changes sign 

- tolerance for root; i.e., if Xi and x 2 are two consecutive approximations for the root, 
the average of x a and x 2 is accepted as a root if 

l*i - x 2 | *s (1 MAX (|x a | MAX |x 2 |)) * <tolerance for root> 

- tolerance for function evaluation where an approximation x is accepted as a root if 
|f(x)| =s <tolerance for the function>. 

Press CONTINUE after each entry. 

4. The root x as well as the functional value f(x) will be printed. 

Example 

User-defined function: 

F = X * LOG(X) - 1 

User entries: 

FIRST GUE:SS«= i 

STEP SIZE= .25 

MAXIMUM NUMBER OF STEPS™ 30 

MAXIMUM NUMBER OF ITERATIONS^ 3D 

TO i... FRANCE FOR R T = :'i. . E •••• 8 

TOLERANCE FOR FUNCTION^ i.E-6 

Results: 

R T === 1.763223 E + F U N C T I ON V A I... U E ~ •••• 3.031967 E •■•■ 8 
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Subprogram Work 

This subprogram is used by subprogram Roota, a root finder subprogram. Work finds the root 
of a user-defined real-valued continuous function by a modified secant method. It requires 
that the functional values of the endpoints of the search interval be of opposite sign. 

Subprogram Utilization 

• File Name - contained in file "Roota", disc 1 

• Calling Syntax 

CALL Work (Root, Err, Tl, Zl, T3, Z3, Itrmxjola, Tolf) 

• Input Parameters 



T1.T3 

Z1.Z3 

Itrmx 

Tola 



Tolf 



Endpoints of search interval. 

f(Tl), f(T3) 

Maximum number of iterations allowed in searching for a 

root. 

Tolerance for the root; i.e., if x : and x 2 are two consecutive 

approximations for the root, the average of x x and x 2 is 

accepted as a root if 

| Xl - x 2 | < (1 MAX (| Xl | MAX |x 2 |)) * Tola. 

Tolerance for the function; i.e., an approximation x is 
accepted as a root if |f(x)| «= Tolf. 



Output Parameters 

Root 
Err 

1 Subprograms Required 

FN F(X) 



Root of f(x) 
f(Root) 



User-defined function which, when given a domain value x, 
returns the value f(x). 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, Zl and Z3 must be of opposite sign. If not, an error 
message is printed and the program pauses. 

"ERROR IN SUBPROGRAM Work." 

"Zl and Z3 ARE NOT OF OPPOSITE SIGN." 

Zl = Z3 = 

You may want to restart the program and correct the data. 

• If subprogram Roota is being used, Work is called only in subprogram Chkit. Chkit speci- 
fically checks that Zl and Z3 are opposite sign. So this error will not occur. 

If subprogram Work is being used independently, then you should check both the logic 
in the calling program and the specific behavior of the function in the interval [Tl, T3]. 

• If Itrmx, the maximum number of iterations is exceeded, an error message is also printed. 

"ERROR IN SUBPROGRAM Work." 

"MAXIMUM # OF ITERATIONS EXCEEDED." 

Kounti = Itrmx - 

The maximum number of iterations may be corrected from the keyboard (e.g., 
Itrmx = <new maximum>, EXECUTE). When CONTINUE is pressed, the program will 
resume, incorporating the new value of Itrmx. There is a reasonable limit here, though. 
Itrmx greater than 50 would probably not make sense. If convergence is not obtained 
after 50 iterations, there is some other problem. You should check to see if f(x) is exhibit- 
ing some unusual behavior in the vicinity of the looked-for root. It might even be neces- 
sary to switch to a simple bisection method, subprogram Bisect, in this subinterval to find 
the root. 

• In general, the modified secant method will be faster and obtain greater accuracy than 
the bisection method. If a large number of roots needs to be located, subprogram Roota 
or Muller should be used. 

Methods and Formulae 

See "Methods and Formulae" section for subprogram Roota on page 1-8. 
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Subprogram Chkit 

This subprogram is used in subprogram Roota. Chkit checks if a root has been found or 
bracketed. 

Subprogram Utilization 

• File Name - contained in file "Roota", disc 1 

• Calling Syntax 

CALL Chkit (Root, Err, Tl, Zl, Step, T2, Z2, Itrmx, Tola, Tolf, Rtn) 

• Input Parameters 



T1.T2 
Z1.Z2 

Step 

Itrmx 
Tola 

Tolf 
Rtn 

Output Parameters 



Endpoints of search interval. 

f(Tl), f(T2). 

Step size for marching when looking for an interval in which 
the function changes sign; step size must be a positive 
value. 

Maximum number of iterations allowed in subprogram 
Work after an interval has been found in which the function 
changes sign. 

Tolerance for the root; i.e., if x 1 and x 2 are two consecutive 
approximations for the root, the average of x Y and x 2 is 
accepted as a root if 

|xj - x 2 | ss (1 MAX (|x : | MAX |x 2 |)) * Tola. 

Tolerance for the function; i.e., an approximation x is 
accepted as a root if |f(x)| =s Tolf. 

Rtn = 1 implies a root was located. 

Rtn = — 1 implies a root was not yet located. 



Root 


Root of f(x 


Err 


f(Root) 


• Subprograms Required 




FN F(X) 


User-define 



Work 



User-defined function which, when given a domain value x, 
returns the value f(x). 

Finds the root of the user-defined real-valued continuous 
function f by a modified secant method. 
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Special Considerations and Programming Hints 

See subprogram Roota (page 1-7 through 1-8). 

Methods and Formulae 

Upon entry into Chkit, the search interval is incremented one step size. A check is performed 
to see if the function value of the new endpoints satisfies the error tolerance, Tolf. If so, the 
root has been located and the subporgram is exited. If not, a check is made to see if the root 
has been bracketed. If so, subprogram Work is called. If not, the subprogram is exited. 



Subprogram Monton 

This subprogram is used in subprogram Roota, a root finder. Monton checks the three points 
(x 1; y x ), (x 2 , y 2 ), (x 3 , y 3 ) to see that y = f(x) is strictly increasing or decreasing. If not, then an 
error message is printed. 

This subprogram insures that the function is well-behaved in the search interval. 

Subprogram Utilization 

• File Name - contained in file "Roota", disc 1 

• Calling Syntax 

CALL Monton (Tolf, XI, X2, X3, Yl, Y2, Y3) 

• Input Parameters 

Tolf Provides tolerance for flatness; also used as tolerance for 

function evaluation in other subprograms. 

XI, X2, X3 Three current values of search interval. 

Y1,Y2,Y3 f(Xl), f(X2), f(X3) 

The input parameters are all left unchanged and there are no output parameters for this 
subprogram. 
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Special Considerations and Programming Hints 

• Upon entering this subprogram, a check is made to see that X1<X2<X3. If not, an error 
message is printed and the program pauses. 

"ERROR IN SUBPROGRAM Monton." 

"ORDINATES XI, X2, X3 ARE NOT IN INCREASING ORDER " 

XI = X2 = X3 = 

If this error occurs, it suggests that the user-defined function is extremely ill-behaved. 
You may want to redefine the function (see page 1-8). 

• Monton' s primary purpose is to check that the function is monotonic in the interval under 
consideration; i.e. x x < x 2 < x 3 or y 1 < y 2 < y 3 . 

There is also cause for alarm if the function is almost flat in the interval; i.e., y l7 y 2 and y 3 
are extremely close in value. If this is the case, then an error message is printed and the 
program pauses. 

"ERROR IN SUBPROGRAM Monton." 

"Y = F(X) IS NOT MONOTONIC OR IS VERY FLAT." 

"SUGGEST TRYING SMALLER STEP SIZE OR DIFFERENT FIRST GUESS " 

XI = X2 = X3 = 

XI = Y2 = Y3 = 

There are many possible causes of this error, the most common being that (X2, Y2) is a 
relative minimum or maximum. Choosing XI or X3 as the first guess may solve this 
problem. Another possibility is to choose a smaller step size. You may want to restart the 
program and make these changes. 

Methods and Formulae 

Monton performs two checks: 

1. that x x < x 2 < x 3 
and 

2. that (y 3 - y 2 ) * (y 2 - y x ) =s Tolf 

This second inequality tests that f(x) is strictly increasing or decreasing in the interval [x 1( x 3 ] 
and that the values y x , y 2 , y 3 are not too close in value. The tolerance for the function, Tolf, is 
used for that purpose here. 
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Subprogram FN Incres 

This subprogram is used in subprogram Roota, a root finder, and is used to check if too many 
marching steps have been taken in searching for a change of sign. 

Subprogram Utilization 

• File Name - contained in file "Roota", disc 1 

• Calling Syntax 

FN Incres (I, Istpmx, XI, X2, X3, Yl, Y2, Y3) 

• Input Parameters 

I Number of steps taken upon entering subprogram. 

Istpmx Maximum number of steps allowed in marching to find an 

interval in which the function changes sign. 

XI, X2, X3 Three points in the interval [x 1; x 3 ] 

Yl, Y2, Y3 Yl - f(Xl), Y2 = f(X2), Y3 = f(X3) 

• Output 

FN Incres I + 1 

Special Considerations and Programming Hints 

• If the maximum number of steps, Istpmx, has been exceeded, an error message is 
printed and the program pauses. 

"ERROR IN SUBPROGRAM Incres." 

"NO CHANGE OF SIGN AFTER <i> STEPS" 

XI = Yl = 

X2 = Y2 = 

X3 = Y3 = 

If you wish to continue with a larger maximum number of steps, Istpmx may be changed 
from the keyboard (e.g., Istpmx = <new maximum>, EXECUTE). Press CONTINUE 
and the program will resume with the new maximum number of steps. If Istpmx = 20 
does not work, you may want to attempt Istpmx = 50. If this is not effective, it would not 
make sense to try Istpmx = 100 or Istpmx = 1000. There is some other difficulty. 
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Muller's Method 

Description 

Given a vector of initial guesses, this subprogram will search for solutions of f(x) = where 
you define the continuous real-valued function f(x). Roots are found by Muller's method, a 
real, quadratic root solver. 

You are required to specify the number of roots to be found, the error tolerance for the 
function evaluations and the number of significant digits desired in the roots, as well as the 
spread criteria for multiple roots and the maximum number of iterations per root. 

Program Usage 

Driver Utilization 

• File Name - "MULLER", disc 1 

The driver "MULLER" sets up the necessary input parameters for the subprogram Muller and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Muller", disc 1 

• Calling Syntax 

CALL Muller (Root(*), Nroots, Itmax, Tolf, Eps, Eta, Digits) 

• Input Paramters 

Root(*)' Vector containing initial guesses of roots. 

Nroots Number of roots to be found. 

Itmax 1 Maximum number of iterations per root. 

Tolf Function tolerance; i.e., x is accepted as a root if 

|f(x)| =s Tolf. 

Eps Spread tolerance for multiple roots. 

Eta Restart value for multiple roots. 



Note 
If the Ith root (Root(I)) has been computed and it is found that 
|Root(I) - Root(J)| =£ Eps where Root(J) is a previously computed 
root, then the computation is restarted with a guess equal to 
Root(I) + Eta. 

Digits Number of significant digits desired in root. 



1 Itmax and Root(») are both input and output parameters. 
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• Output Parameters 

Root(*) 1 Vector containing roots of the function. 

Itmax 1 Number of iterations required to find final root. 

• Subprograms Required 

FN F(X) User-defined function which, when given a domain value x 

returns the value f(x). 



Special Considerations and Programming Hints 

• Upon entry into subprogram Muller, there is a bad data check. If the check detects any 
"nonsense" data, an error message will be printed and the program pauses. 

"ERROR IN SUBPROGRAM Muller." 
Nroots = Itmax = 

Tolf = Digits = 

The data may be corrected from the keyboard (e.g., Tolf = IE - 8, EXECUTE). When 
CONTINUE is pressed, the program will resume execution at the next line. 

• If the maximum number of iterations, Itmax, is exceeded in searching for any one root, I, 
an error message will be printed. 

"ERROR IN SUBPROGRAM Muller." 

"MAXIMUM # OF ITERATIONS EXCEEDED ON Root <i>" 

Root(I) is assigned the value 9.999999E99. The subprogram continues execution with a 
search for the next root. 

• To define the function f(x): 

1. PAUSE the current program (if one is running). 

2. Type SCRATCH and press EXECUTE. 

3. In EDIT mode, type in the following function subprogram pressing ENTER after 
each line: 

10 DEF FNF(X) 

20 F: F = < user-defined function of the form (X-3)*(X + 4)> 

30 RETURN F 

40 FNEND 

4. RE-STORE the function subprogram on file "KRYWEN". 

5. LOAD file "MULLER" and press RUN. 

(See "Redefining User Functions" in the Introduction, following the Table of Contents, for 
another way to modify the function f(x). ) 
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• Muller's technique has the quickest rate of convergence of all the methods contained in 
this section. If you have good initial estimates of the roots, this method is the one to use. 

However, there are times this method gives misleading results. For example, if the first 
fifty positive roots of a particular function are desired and you do not have good initial 
estimates, Muller may not generate the desired results. Some roots may be skipped over 
or some negative roots may be generated. 

For example, assume we want to locate the first 6 positive roots of f = sin(x) and we 
have poor initial estimates: 

Example 

User-defined function: 

F = SIN(X) 



User entries: 



11: of root!;; to be found^ s 

maximum * of" iterations- 20 

tolerance foe function-^ i . £■ -6 

op read tolerance for multiple roots'-- :"i. . 

rest art v'alue for multiple roots ^ .0001 

* of' significant digits in r00t= 6 



INITIAL GUESS i i) 

INITIAL GUESS ( 2) 

INITIAL GUESS ( 3) 

INITIAL GUESS ( 4"r- 

INITIAL GUESS ( S) 

INITIAL GUESS ( 6): 



OOOOOOEv-00 
r *■ 
0O0OO0E+OO 
E •!- (.1 
00 00 (IE +00 
OOOOOOE*- 



Results: 



ROOTS 



R (j ci t ( 
R o o t ( 
R o o t ( 
R (i o t ( 
E oot( 
R o o !' ( 



i)-"- 2.620 439E 08 
2)---~ 3. 14lS93Ev-0 
3 ) ■■-■ i= On 51 BSE MM) 
A j ■■■■■■■■ 9. 4 24 7 78 E + 
3)= 1 . 25S637EMH 
&)-~-3.i4i5?3EM)0 
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On the other hand, assume that our initial guesses are good: 
User-defined function: 

F = SIN(X) 



User entries: 



Results: 



# OF ROOTS TO BE FOUND-- 6 
MAXIMUM t OF ITERATIONS- 20 

T O L E R A M C E F R I"' U N C T 1 N :::: i . E -6 

SPREAD TOLERANCE FOR MULTIPLE ROOTS- i . E-8 

R E S T ART V A I... U E FOR h U L T I P L E R T S :::: . 1 

# OF SI ON IF 10 ANT DIGITS IN ROOT= b 



INITIAL. GUESS ( 1. > ~ 3 . 00 OEM) 

I N I T I A L G U E S S < 2. ) - 6.0000 E + 00 

I N I T I A L G U E S S ( 3 ) = 9.00 E + 

I N I T I A I... G U E 8 S ( 4 ) -= :1. . 2 E + 1 

I N I T I A !... QUE S S < S) :::: 1 . 5 E + 1 

I N I T I A I... G U E S S < 6 ) ■■- i . 8 E + 1 



ROOTS 



Rod I ( i ) -" 3 . 5 4 :i. S 9 3 E + 

RootC 2)-== 6.283185E+0 

R oot( 3>™ 9 . 424778E+0 

Roo'K 4)= 1 . 2366 37E+0 i 

Root ( 3) "" 1 . 370796E+01 

Root( 6) :: " 1 .8 84 9 86 EMU 
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Methods and Formulae 

Given a vector of initial guesses, subprogram Muller will search for solutions of f(x) = 
where you define the continuous real-valued function f(x). 

The algorithm is based on an extension of Muller's method by Werner L. Frank. This proce- 
dure does not depend on any prior knowledge of the location of the roots nor on any special 
starting process. All that is required is the ability to evaluate f(x) for any desired value of x. 
Multiple roots can also be obtained. 

Given that you supply an initial guess, Root(I), 
P = .9 * Root(l) 
Pi = 1.1 *Root(I) 
P 2 = Root(I) if Root(I) 40. 

If Root(I) = 0, P = - 1, P x = 1, and P 2 = 0. From these three starting values, you choose 
Rt, the next approximation to the root, as one of the zeros of the second degree polynomial 
which passes through the functional values f(P), f(P x ) and f(P 2 ). Successive approximations 
are obtained by repeating the quadratic fit: 



where d if3 = 



Rt, + 3 = Rt i+ 2 + (Rt i+ 2 - Rt i + 1 ) * d, + 3 
- 2*F(Rt i + 2 )*(l + d i + 2 ) 



b i + 2 ± b 2 i + 2 - 4* F(Rt i + 2 ) * d i + 2 * (1 + d i + 2 )*{ F(Rti) * d i + 2 - F(Rt i + 1 ) * (1 + d i + 2 ) + F(Rt i + 2 )} 1/2 

andbi + 2 = F(Rt,)*d 2 +2 - F(Rt i + 1 )*(l + d j + 2 ) 2 + F(Rt i + 2 ) * (1 + 2*d i + 2 ) 

The sign in the denominator is chosen to give d i + 3 the smaller magnitude. 

Having found (R - 1) roots, you determine the Rth root by solving the equation f R (x) = 
where: 

f R (x) = ^ (R = 2,3,...). 

R - 1 

n (x - x,) 
i = 1 

The process halts when successive iterants pass one of the two convergence tests: 

|H| < |Rt| *1(T Di3its where H - Rt - previous Rt 

Digits — number of significant digits 
desired in root 
or 

|f(Rt)| < Tolf and |f(Rt)| < Tolf where Tolf is the desired function tolerance. 

If the maximum number of iterations is exceeded on any root, an error message will be 
printed and the root will default to the value 9.999999E99. 
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References 

1. D.E. Muller, "A Method for Solving Algebraic Equations Using an Automatic Compu- 
ter", MTAC 10 (1956). 

2. W.L. Frank, "Finding Zeros of Arbitrary Functions", J.ACM 5 (1958). 

3. B. Leavenworth, "Algorithm 25: Real Zeros of an Arbitrary Function", Comm. of ACM 
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User Instructions 

1. If you have not already defined the function you wish to solve, you should do so at this 
time. 

2. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "MULLER" and press RUN. 

3. You will be asked to supply entries for the following items: 

- number of roots to be found 

- maximum number of iterations permitted per root 

- tolerance for functional evaluation where x is accepted as a root of |f(x)| =s <this 
tolerance > 

- spread tolerance for multiple roots. If the Ith root, Root(I), has been computed and it 
is found that |Root(I) - Root(J)| =£ <this tolerance> where Root(J) is a 
previously computed root, then the computation is restarted with a guess equal 

to Root(I) + <restart value>. 

- restart value for multiple roots 

- number of significant digits in root for it to be acceptable 

- initial guesses for the number of roots expected. 
Press CONTINUE after each entry. 

4. Values will be printed for all roots found. The value 9.999999E99 is inserted for any 
roots not found. 



1-22 Root Finders 



Example 

User-defined function: 

F = SIN(X) 
User entries: 



t Of ROOTS TO BE: I'" OUHI):::: 4 

nAXIMUM * OF ITERATIONS^ 20 

TO i... F R A N E F P F 1M "f 1 N ~ i . E ■■■■ 8 

SPREAD TOL.FRAMCF FOR MULTIPLE ROOT'S^ i. F ■••• i. 

RESTART UAU..1F FOR HL1LTIPLF ROOTS'™ .00-01 

* OF SIGNIFICANT DIGITS IN R!'JOT~ 8 



I N I T I A I... G U E S S >; :i ) ■■■■ . () E •*• 

i n i r :i: al guess < 2 :> ■■■■■■■ ;■» c e-<- 

I. N I. T I A!... G OESS ( 3 ) "- 4.00 E -:- 

i N I T I A 1... G I! E S S < 4 ) ■■■■■ 6 . s: + 



Results: 



:GOT 



R o o -i; ( i ) ■"■ . £ f () 

Roof ( 2) ■■■■■■■■ 3. :'i 41593EX)!) 

R oof ( 3 ) ~ t , 283:1 BSE •*•() 

Roo t ( 4) - 9 . 4247781:; •* 
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Siljak's Method 

Description 

This subprogram will find all roots, Z, of polynomials of the form 

a + ib + (a x + ibJZ + (a 2 + ib 2 )Z 2 + — + (a n + ib n )Z n = 0. 

Roots are found by expressing the polynomials in terms of Siljak functions and using the 
method of steepest descent to determine the zeros. 

You are required to enter the complex coefficients of the polynomial as well as its degree. 
Tolerances for the root and for function evaluations and the maximum number of iterations 
also need to be entered. 

Program Usage 

Driver Utilization 

• File Name - "SILJAK", disc 1 

The driver "SILJAK" sets up the necessary input parameters for the subprogram Siljak and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name -"Siljak", disci 

• Calling Syntax 

CALL Siljak (N, Rcoef(*), Icoef(*), Tola, Tolf, Itmax, Rroot(*), Iroot(*)) 

• Input Parameters 

N Degree of polynomial; number of roots to be found. 

Rcoef(*) Vector containing the real parts of the coefficients of the 

polynomial where the subscript corresponds to the expo- 
nent of the variable; subscripted from to N; e.g., for 
a + ib +(a x + ib x )Z + (a 2 + ib 2 )Z 2 + --- 
+ (a n + ib n )Z n = 0, Rcoef(3) would be the real part of the 
coefficient of the Z 3 term, a 3 . 

Icoef(*) Vector containing the imaginary parts of the coefficients of 

the polynomial where the subscript corresponds to the ex- 
ponent of the variable; subscripted from to N. 

Tola Tolerance for the root. 

Tolf Tolerance for the function evaluation. 

Itmax Maximum number of iterations permitted in searching for 

any one root. 

• Output Parameters 

Rroot(*) Vector containing the real parts of the roots of the polyno- 

mial; subscripted from 1 to N. 

Iroot(*) Vector containing the imaginary parts of the roots of the 

polynomial; subscripted from 1 to N. 
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Special Considerations and Programming Hints 

• Upon entry into subprogram Siljak, there is a bad data check. If the check detects "non- 
sense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Siljak." 
N = Tola = 

Tolf = Itmax = 

The data may be corrected from the keyboard (e.g., Tola = IE - 6, EXECUTE). When 
CONTINUE is pressed, the program will resume execution at the next line. 

• The subprogram has been set to quarter the interval size at most 20 times. If this max- 
imum is exceeded, an error message is printed: 

"ERROR IN SUBPROGRAM Siljak." 

"THE INTERVAL SIZE HAS BEEN QUARTERED 20 TIMES AND THE TOLER- 
ANCE FOR FUNCTIONAL EVALUATIONS IS STILL NOT MET." 
Tolf = U = V = 

When CONTINUE is pressed the subprogram is exited with all known information. All 
roots not found will contain the value 9.999999E99. 

If more than 20 quarterings are required, there is probably some unusual behavior in the 
function. 

• If the maximum number of iterations has been exceeded, the following error message 
will be printed and Siljak will pause: 

"ERROR IN SUBPROGRAM Siljak." 

"MAXIMUM # OF ITERATIONS HAS BEEN EXCEEDED." 

L = Itmax = 

When CONTINUE is pressed, the subprogram will be exited. Again, all roots not located 
will contain the value 9.999999E99. 

There are two possible solutions here. First, increase Itmax, the maximum number of 
iterations allowed; or second, allow larger error tolerances, Tola and Tolf. The iteration 
counter, L, is reset after each root is located, so attempting to solve a polynomial of large 
degree should not cause this error. 
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Methods and Formulae 

Roots are found by expressing the polynomial in terms of Siljak functions and using the 
method of steepest descent to determine the zeros. Once a root is found, the polynomial is 
reduced by synthetic division and the process is repeated. The last root is computed alge- 
braically. The algorithm is very accurate and stable; it will virtually always find the roots and 
you are not required to provide an initial value. Multiple roots are found at some slightly 
reduced accuracy, and higher order polynomials may show some loss of accuracy as more 
roots are found. In general, the program will find "normally" spaced roots accurate to better 
than 6 decimal places. Newton's method could find the roots faster but convergence is not 
guaranteed and with Siljak's method, no a priori information such as the derivative is neces- 
sary. 

f(Z) = 2 (a k + ib k )Z k = 

k = 

Siljak functions x k and y k are defined by Z k = x k + iy k and may be calculated recursively 
where x + iy are the root approximations. 



x = 1, x x = .1, y = 0, yi = 1 
x k + 2 = 2xx k + 1 - (x 2 + y 2 )x k 
y k + 2 = 2xy k + 1 - (x 2 + y 2 )y k 



M- = 2 (a k x k - b k y k ) 

k = 



d|x 

dx - ^ Kia k A k _! - u k y k _ x / 
k= 1 



2 k(a k x k _i- b k y k _ ; 



3x 



== 2 k(a k y k _a+ b k x k _j) 
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Reference 

1. Moore, J.B., "A Convergent Algorithm for Solving Polynomial Equations", Journal of 
the Association for Computing Machinery, Vol. 14, No. 2 (April, 1967), pp. 311-315. 

User Instructions 

1. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "SILJAK" and press RUN. 

2. You will be asked to supply entries for the following items: 

- degree of polynomial (this is also the number of roots that will be found) 

- maximum number of iterations allowed per root 

- tolerance for roots where a root is accepted if the difference in value of the root 
approximations of two successive iterations is less than this tolerance 

- tolerance for functional evaluations where the root x is accepted if |f(x)| =£ <this 
tolerance> 

- the real part and the imaginary part for each coefficient of the polynomial (the array 
subscript shown corresponds to the exponent of the variable). 

Press CONTINUE after each entry. 

3. The real and imaginary parts of the roots will be printed. Any roots not found will 
contain the value 9.999999E99. 



Example 

User entries: 



... Q 



DEGREE OF POLYNOMIAL™ 4 
MAX # OF ITERATIONS™ 21) 

tolerance: for roots- i . i 

TOLERANCE FOR FUNCTIONAL EVALUATIONS™ i 



C E F F I C I E N T S : I! ( R c o e f < ) + 1 <:: o e f ( ) * I ) + 

(Rcoefd ) + koef ( i)#'.i:)*X A i + . .1 

REAL IMAGINARY 

-■ 1 . 60 0E + 1 0.0 QE + 

. 00OOO0E+00 . 00OO0OE+O0 

. 0OOO00E+O0 . OOOOOQE+0 

. 00O000E+O0 . 00 00 0E+O0 

1 E ■'<■ . E+ 



Results: 



ROOTS 



EAL IMAGINARY 



•2 . OOOOOOE+OO -1 . 747497E-24 

2. 000 000E+0 :i. . 747497E -24 

i.74 8 ;:> 4 E - 2 4 •••■ 2.00 E + 

i. .748 65E-24 2. OOOOOOE + 
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Chapter A 
Integration 



Introduction 

Description 

This section contains routines for numerically evaluating integrals. There are many possible 
reasons for using numerical integration, three of which are: 

1. the analytical form of a simple integral can be quite complicated. For example, 



J dx / V (x 2 (x 2 + 25) a ) = - 1 / (50x 2 V (x 2 + 25)) - 3 / (2(5 4 )Vx 2 + 25) + 
(3 / (2(5 5 ))) log (5 + Vx 2 + 25/x) 

The evaluation of this expression may involve more calculations than evaluation by 
numerical methods. 

_ 2 

2. many integrals cannot be expressed in finite form. For example, j e x dx. 

3. the integrand may not be known explicitly, but may be expressed as a set of collocation 
points. 



Routines 

• Simp - the easiest method to use for well-behaved functions. 

• Filon - for functions of the form f(x)sin(x) or f(x)cos(x). 

• Cadre - used for a great many functions when other methods fail. 

• Inteq - numerical integration with equally-spaced data points. 
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Simpson's Rule 

Description 

This subprogram approximates | f(x)dx for the user-defined continuous function f(x) on the 
interval [a,b]. You are required to specify the maximum number of iterations permitted and 
the error tolerance between successive calculations of the integral. 

Program Usage 

Driver Utilization 

• File Name -"SIMP", disci 

The driver "SIMP" sets up the necessary input parameters for the subprogram Simp and 
prints out the resulting value of the integral. Intermediate integral values may also be printed if 
you wish. 

Subprogram Utilization 

• File Name - "Simp", disc 1 

• Calling Syntax 

CALL Simp (Low, Up, Int, Fig, Itmax, Tol) 

• Input Parameters 

Low Lower bound of interval of integration. 

Up Upper bound of interval of integration. 

Fig Should intermediate results be printed? 

Fig = 1 implies yes 
Fig = - 1 implies no 

Itmax Maximum number of interval halvings. 

Tol Error tolerance between successive calculations of integral. 

• Output Parameters 

Int Value of integral. 

• Subprograms Required 

FN F(X) Function to be integrated must be defined in function sub- 

program FN F(X). 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, a bad data check is performed. If the program detects 
"nonsense" data, the program will pause and the following error message will be 
printed: 

"ERROR IN SUBPROGRAM Simp." 
Up = Low = Fig = 

Itrnax = Tol = 

The data may be corrected from the keyboard (e.g., Tol = IE - 4, EXECUTE). When 
CONTINUE is pressed, the program will resume execution at the next line. 

• If the maximum number of iterations is exceeded, the following error message will be 
printed: 

"ERROR IN SUBPROGRAM Simp." 
"MAXIMUM # OF ITERATIONS EXCEEDED." 
Int = Intold = 

Tol = Itmax = 

At this point, the tolerance or the maximum number of iterations may be increased from 
the keyboard (e.g., Itmax = 30, EXECUTE) or the value of the integral may be accepted 
as is. If the data is corrected from the keyboard, pressing CONTINUE will cause execu- 
tion to resume. 



• To define the function f(x): 

1. PAUSE the current program (if one is running). 

2. Type SCRATCH and press EXECUTE. 

3. In EDIT mode, type in the following function subprogram pressing ENTER after 
each line. 

10 DEF FNF(X) 

20 F: F = <user-defined function of the form (X-3)*(X + 4)> 

30 RETURN F 

40 FNEND 

4. RE-STORE the function subprogram on file "KRYWEN". 

5. LOAD file "SIMP" and press RUN. 

(See "Redefining User Functions" in the Introduction, following the Table of Contents, for 
another way to modify the function f(x).) 
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Methods and Formulae 

This subprogram will approximate f b f(x)dx for the user-defined function f(x) which is defined 
in a function subprogram. The function must be continuous over the interval [a,b]. 

Simpson's one-third rule: 

,» h 

I f(x)dx - - [f(a) + 4f(a + h) + 2f(a + 2h) + 4f(a + 3h) + - + 4f(a + (n - l)h) + f(a + nh)} 

where n = number of intervals, 
h = - 1 — zr^- = interval size 



Reference 

1. Beckett, Royce and Hurt, James, Numerical Calculations and Algorithms , New York: 
McGrawHill, 1967, pp. 166-169. 

User Instructions 

1. If you have not already defined the function you wish to solve, you should do so at this 
time. 

2. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "SIMP" and press RUN. 

3. You will be asked to supply entries for the following items: 

- lower bound of the interval 

- upper bound of the interval 

- maximum number of interval halvings 

- error tolerance where the value of the integral is accepted if the difference 
in value of two successive approximations is less than this tolerance. 

Press CONTINUE after each entry. 

4. The program will print the value of the integral as well as intermediate data points if 
they were requested. 
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Example 

User-defined function: 

F - X*X*SIN(3*X) 

User entries: , ...,,,...,., ,..,.,,, , 

U i :: ' P E R BOUND = 1 . 4 7 i 9 7 5 5 1 2 

MAX # Or INTERVAL HALVINGS^ 10 

ERROR TOLERANCE* . 0.1 

Intermediate data points: .,. . ,... ... , , . . .... ,., .... .., , . .. , r , , ... ,., ,, , •. ,••■ n - 

i- :i: n t i;:: r v a i... -- 4 i n t e g r a i... -■ 2 . ± ? 2 1 6 1;;: ■•■• .1. 

# INTERNAL"- 8 'INTEGRAL-- 2 . 17379SE-' i 

* INTERVAL- 16 INTEGRAL. = 2 . 1 7392 iE-0 i 

Results: ... ., ,.. , , 

IN 1 H/..RH !...- 
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Filon's Method 

Description 

This subprogram approximates |f(x)sin(tx)dx or f f(x)cos(tx)dx for the user-defined con- 
tinuous function f(x) on the interval [a,b]. Both integrals are found using Filon's method, a 
special case of Simpson's rule for oscillating functions. 

You are required to specify an array of function values, the number of points of evaluation, 
the number of oscillations (t), and the endpoints of the interval. 

Program Usage 

Driver Utilization 

• File Name - "FILON", disci 

The driver "FILON" sets up the necessary input parameters for the subprogram Filon and 
prints out the resulting value of the integral. 

Subprogram Utilization 

• File Name - "Filon", disc 1 

• Calling Syntax 

CALL Filon (F(*), T, A, B, Nb, Key, Int) 

• Input Parameters 



F(*) 
T 

A, B 

Nb 

Key 

• Output Parameter 

Int 

• Subprograms Required 

Cnsts (Alpha, 
Beta, Gamma, Theta) 



Vector containing table of functional values calculated at Nb 
equidistant points from A to B; subscripted from 1 to Nb. 

Number of oscillations of the trigonometric function; i.e., 
sin(tx) or cos(tx). 

Endpoints of interval. 

Number of equidistant points from A to B; Nb must be odd 
and Nb > T + 1. 

Key = 1 implies cosine integral is evaluated. 
Key = - 1 implies sine integral is evaluated. 



Value of the integral. 



Subprogram which computes the three constants for use in 
the Filon integration formula. 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, a bad data check is performed. If the program detects 
"nonsense" data, the program execution will pause and the following error message will 
be printed: 

"ERROR IN SUBPROGRAM Filon." 
T = A = B = 

Nb = Key = 

The data may be corrected from the keyboard (e.g., Nb = 71, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• In order to use Filon' s formula, an odd number of points, Nb, is required. If the subpro- 
gram detects an even number of points, an error message is printed and execution 
pauses. 

"ERROR IN SUBPROGRAM Filon." 
"ODD # OF POINTS REQUIRED." 

Nb = 

Again, the data may be corrected from the keyboard. When CONTINUE is pressed, the 
program will resume execution at the next line. 

• For a good approximation of the integral, Theta values (Theta = t*h) should be kept less 
than 1. So, for example, if you want to integrate |x 2 sin(10x)dx, t = 10. Nb should be 
chosen greater than 11 since: 



lOh < 1, h < 



(b ~ a) < J I < J — , Nb > 11. 



10 Nb - 1 10 Nb - 1 10 



Care should also be taken not to choose too large an Nb since this increases the running 
time of the program. 

• The vector F(*) must be dimensioned in the calling program: DIM F(l:Nb) where Nb is 
the number of equidistant points of the interval of integration. 

• To define the function f(x) where f(x)sin tx or f(x)cos tx is to be integrated: 

1. PAUSE the current program (if one is running). 

2. Type SCRATCH and press EXECUTE. 

3. In EDIT mode, type in the following function subprogram pressing ENTER after 
each line. 

10 DEF FNF(X) 

20 F: F = <user-defined function of the form (X - 3) * (X + 4) > 

30 RETURN F 

40 FNEND 

4. RE-STORE the function subprogram on file "KRYWEN". 

5. LOAD file "FILON" and press RUN. 

(See "Redefining User Functions" in the Introduction, following the Table of Contents, for 
another way to modify the function f(x).) 
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Methods and Formulae 

This subprogram computes | f(x)cos(tx)dx or { b f(x)sin(tx)dx. Both integrals are found using 
Filon's method, a special case of Simpson's rule. 

It is assumed that F(*) is an array of functional values calculated at Nb (odd) equidistant 
points from a to b and h is the interval size, 

h = A^L 

Nb - 1 
| b f(x)sin(tx)dx - h[aS a + (3S P + yS y ] 

| b f(x)cos(tx)dx - h[aC a + pCp + yC y ] 

where a= [t 2 h 2 + th sin(th)cos(th)- 2 sin(th) 2 ] / (th) 3 
3= 2[th(l + cos(th) 2 - 2 sin(th)cos(th)] / (th) 3 
7= 4[sin(th) - th cos(th)] / (th) 3 
S a = - [f(b)cos(bt) - f(a)cos(at)] 

Nb-l 



2 



S 3 = .5[f(a)sin(at) - f(b)sin(bt)] + 1 f(2, + l)sin[(a + 2ih)t] 



i=i 

Nb-l 



2 

S 7 = 2 f(2i)sin[(a + (2i - l)h)t] 

i=i 

C a = f(b)sin(bt) - f(a)sin(at) N b-i 

2 

C p = .5[f(a)cos(at)- f(b)cos(bt)] + X f(2i + l)cos[(a + 2ih)t] 

i=i 

Nb-l 



2 

C 7 = 2 f(2i)cos[(a + (2i - l)h)t] 



i = i 



User Instructions 

1. If you have not already defined the function you wish to solve, you should do so at this 
time. 

2. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "FILON" and press RUN. 

3. You will be asked to supply entries for the following items: 

- number of points of integration (odd number required) 

- form of integrand, whether f(x)sin(tx) or f(x)cos(tx) 

- lower bound of the interval 

- upper bound of the interval 

- frequency or number of oscillations of the trigonometric function 
(must be greater than 0) 

Press CONTINUE after each entry. 

4. The program will print the value of the integral. 
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Examples 

User-defined function: 

F = X*X 
User entries: 



* OF POINTS OF INTEGRATION- 31 

Key- i 

LOWER BOUND= 

UPPER BOUND™ 1.5707963268 

FREQUENCY" i 



Result: INTEGRAL- 4 . 67401 1E-01 



User-defined function: 

F = LOG(l + X) 

User entries: i( 0| ... p0INTS Qf INTEGRATION- 3i 

l< e y -■■ - i 

LOWER BOUND™ 

UPPER BOUND" 6.2831853072 

FREQUENCY™ 10 

6SU t: IN T E G R A L •■= - 1 . 9 7 6 6 4 E - i 

Using the same function (F = LOG(l + X)): 

User entries: ( , Qf pQ; , NT8 0f , INTEGRATION- 71 

Key - •••■ i 

LOWER BOUND- 

UPPER BOUND- 6. 283 1853072 

FREQUENCY™ 10 



Results: 

:i:NTEGRAL™--l.?76264E"-0i 
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Cautious Adaptive Romberg Extrapolation 

Description 

This subprogram uses cautious adaptive Romberg extrapolation to approximate / f(x)dx for 
the user-defined function f(x) on the interval [a,b]. 

You are required to specify the lower and upper bounds, the relative and absolute error 
tolerances and the user-defined function subprogram. 

Program Usage 

Driver Utilization 

• File Name - "CADRE", disc 1 

The driver "CADRE" sets up the necessary input parameters for the subprogram Cadre and 
prints out the resulting value of the integral, the estimated error term, and the suggested 
reliability of the results. 

Subprogram Utilization 

• File Name - "Cadre", disc 1 

• Calling Syntax 

CALL Cadre (A, B, Aerr, Rerr, Err, Fig, Cadre) 

• Input Parameters 

A Lower bound of interval of integration. 

B Upper bound of interval of integration. 

Aerr Desired absolute error in the answer. 

Rerr Desired relative error in the answer; Rerr must be in the 

range (0,0.1); e.g., Rerr = 0.1 indicates that the estimate of 
the integral is to be correct to one digit, whereas 
Rerr = IE - 4 calls for four digits of accuracy. 



Integration 2-11 



• Output Parameters 

Err 
Fig 



Cadre 
Subprograms Required 
FN F(X) 



Estimated bound on the absolute error of the integral. 
An integer between 1 and 5 indicating what difficulties were 
met with specifically. 
Fig = 1, all is well 

Fig = 2, one or more singularities were successfully 
handled 

Fig = 3, in some subintervals, the estimate Vint was 
accepted merely because Err was small even though no reg- 
ular behavior could be recognized. 
Fig = 4, failure, overflow of stack Ts. 

Fig = 5, failure, too small a subinterval is required. This 
may be due to too much noise in the function (relative to 
the given error requirements) or due to a poorly behaved 
integrand. 
Value of the integral. 



Function to be integrated must be defined in the function 
subprogram FN F(X). 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, a bad data check is performed. If the subprogram 
detects "nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Cadre." 

A = B = Aerr = Rerr = 

The data may be corrected from the keyboard (e.g., Aerr = IE - 6, EXECUTE). When 
CONTINUE is pressed, the program will resume execution at the next line. 

• Cadre can, in many cases, handle jump discontinuities and certain algebraic discon- 
tinuities. See reference 1 for full details. 

• The relative error parameter Rerr must be in the interval [0,0.1]. For example, 
Rerr = 0. 1 indicates that the estimate of the integral is to be correct to one digit, whereas 
Rerr = IE - 4 calls for four digits of accuracy. 

• The absolute error parameter, Aerr, should be nonnegative. In order to give a reasonable 
value for Aerr, you must know the approximate magnitude of the integral being com- 
puted. In many cases, it is satisfactory to use Aerr = 0. In this case, only the relative 
error requirement is satisfied in the computation. 

• DeBoor's original program has been modified in a number of places. Most error mes- 
sages have been deleted and the size of some of the arrays was shortened. Reference 1 
contains the original version of the algorithm. 

• The variable Fig indicates the suggested reliability of the solution. To quote DeBoor: 

"A very cautious man would accept Int only if Fig is 1 or 2. The 
merely reasonable man would keep the faith even if Fig is 3. The 
adventurous man is quite often right in accepting Int even if Fig is 4 
or 5." 

• To define the function f(x): 

1. PAUSE the current program (if one is running). 

2. Type SCRATCH and press EXECUTE. 

3. In EDIT mode, type in the following function subprogram pressing ENTER after 
each line: 

10 DEF FNF(X) 

20 F: F = <user-defined function of the form (X-3) * (X + 4)> 

30 RETURN F 

40 FNEND 

4. RE-STORE the function on subprogram on file"KRYWEN". 

5. LOAD file "CADRE" and press RUN. 

(See "Redefining User Functions" in the Introduction, following the Table of Contents, 
for another way to modify the function f(x).) 
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Methods and Formulae 

The subprogram Cadre attempts to solve the following problem: Given the real function f, two 
real numbers a and b, and two nonnegative numbers Aerr and Rerr, find a number Int such 
that 

I ; b f(x)dx - Int | =s = MAX(Aerr, Rerr * | f b f(x)dx |). 

For this, the subprogram employs an adaptive scheme, whereby Int is found as the sum of 
estimates for the integral of f(x) over suitably small subintervals of the given interval of in- 
tegration. Starting with the interval of integration itself as the first such subinterval, Cadre 
attempts to find an acceptable estimate on a given subinterval by cautious Romberg extrap- 
olation. If this attempt fails, the subinterval is divided into two subintervals of equal length, 
each of which is now considered separately. For the sake of economy, values of f(x), once 
calculated, are saved until they are successfully used in estimating the integral over some 
subinterval to which they belong. For a more detailed description of this algorithm, see refer- 
ence 1. 

Reference 

1. DeBoor, Carl, "CADRE: An Algorithm for Numerical Quadrature". Mathematical Soft- 
ware (John R. Rice, Ed.), New York: Academic Press, 1971, Chapter 7. 

User Instruction 

1. If you have not already defined the function you wish to solve, you should do so at this 
time. 

2. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "CADRE" and press RUN. 

3. You will be asked to supply entries for the following items: 

- lower bound of the interval of integration 

- upper bound of the interval of integration 

- absolute error tolerance 

- relative error tolerance. 

Press CONTINUE after each entry. 

4. The program will print the integral value, the estimated error term, and the suggested 
reliability of the results. 
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Examples 

User-defined function: 

F = X*X*SIN(3*X) 
User entries: 



.OWER BOUND- 
li i :: ' P E R B G i J M D -" :i. . A 7 A 9 7 5 '"> i Z 
ABSOLUTE ERROR™ A.P-iO 
RELATIVE ERROR- 1 . E--8 

Results: ,.. , 

J.N SEGRAi...™ 2 1 7 392 BE -01 

E S T .1' M A T E I) E R R R = 4.55782;? E - 1 i 
FL.G-i 

User-defined function: 

F = LOG(EXP(l)/X) 

User entries: 

LOWER BOUND- 1 . E~:U 

UPPER POUND™ :'i. 
ABSOLUTE ERRORS i . E-S ■ 



Results: 

INTEGRAL- 2. Si 0E + 

E B T I ri A T E D E R R R -" 2.419 8 S 7 E - 6 
PEG- 2 
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Numerical Integration 
with Equally-Spaced Data Points 

Description 

This subprogram approximates / f(x)dx where f(x) is represented by discrete function values 
f(Xj) at equally-spaced points x v 

You are required to supply the increment between intervals, the number of data points (which 
must be odd), and the functional value at each of the data points. 

Program Usage 

Driver Utilization 

• File Name - "INTEQ", disc 1 

The driver "INTEQ" sets up the necessary input parameters for the subprogram Inteq and 
prints out the resulting value of the integral. 

Subprogram Utilization 

• File Name - "Inteq", disc 1 

• Calling Syntax 

CALL Inteq (N, Inc, F(*), Int) 

• Input Paramters 

N Number of data points; this must be odd. 

Inc The increment between equally-spaced data points. 

F(*) Vector containing the functional values in increasing order 

of domain value; subscripted from 1 to N. 

• Output Parameters 

Int The integral over the interval [F(l), F(N)]. 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, there are two error checks. If there are not at least 
three data points, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Inteq." 

"AT LEAST THREE DATA POINTS ARE REQUIRED." 

N = 

If there is not an odd number of data points, the following error message is printed and 
the program pauses: 

"ERROR IN SUBPROGRAM Inteq." 

"ODD # OF DATA POINTS IS REQUIRED." 

N = 

The number of data points can be corrected from the keyboard (e.g., N = <new num- 
ber of data points>, EXECUTE). When CONTINUE is pressed, the program will con- 
tinue execution at the next line. 

• Vector F(*) must be dimensioned in the calling program: DIM F(1:N) where N is the 
number of data points being entered. 

Methods and Formulae 

This subprogram approximates /f(x)dx where discrete values of f(x) are known at equally- 
spaced points over the interval [a,b]. Simpson's one-third rule is used to approximate the 
integral. 

| b f(x)dx * -{f(a) + 4f(a+ h) + 2f(a + 2h) + 4f(a + 3h)+ ... + 4f(a + (n- l)h) + f(a + nh)} 
where: n = number of intervals (number of data points minus one) 



Reference 

1. Beckett, Royce and Hunt, James, Numerical Calculations and Algorithms, New York: 
McGraw-Hill, 1967, pp. 166-169. 



User Instructions 

1. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "INTEQ" and press RUN. 

2. You will be asked to supply entries to the following items: 

- number of data points (odd number required) 

- increment between data points 

- appropriate function value of each data point. 
Press CONTINUE after each entry. 

3. The program will print the value of the integral. 
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Examples 

User entries: # Q| , DAJA po;|:NTS ,, i± 
INCREMENTS . i 

FUNCTION VALUES-. 



Result: 



F < i ) :::: 


. 


0O0OOOE 


+ 


F< 2> :::: 


i . 


OOF. 


-Oi 


F( 3) = 


2. 


00 0E 


-Oi 


F ( 4 ) = 


3 


0OO0OOE 


-Oi 


F( 5_) = 


4. 


0O0000E 


•■•■ 1 


F ( h ) ~ 


5 


0O0OO0E 


-Oi 


F( 7) :::: 


6 


. OOOOOOE 


•••• i 


F( 8) = 


7 


. 0OOOO0E 


-Oi 


F( 9)== 


8 


. OO00O0E 


-Oi 


f',111) - 


9 


. 0O00O0E 


-Oi 


F ( i i ) = 


1 


. OOOOOOE 


+ 


INTEGRAL" 


5 


, OOOOOOE 


■••■ 1 



User entries: ,....,..,...,.,.., 

* UF DA) A Pl!INIb= 

INCREMENT- .23 



FUNCTION VALUES • 



Result: 



. 00 00 OE+00 
2.80 00 OE + 
3. 80 00 OE+0 

s . ;?o ooo OE+oo 

7 E + 
9 . cO00O0E + 
1 . 2100 0E+Oi 
i . 56000 OE+Oi 
2 . E + i 



INTEGRAL- i . 6 4.1.667E + i. 
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Numerical Integration 
with Unequally-Spaced Data Points 

Description 

This subprogram approximates / f(x)dx where f(x) is represented by discrete functional values 
for unequally-spaced domain values x over the interval |a,b]. You are required to input the 
data points and the tolerance desired. 

Program Usage 

Driver Utilization 

• File Name- "INTUN", disc 1 

The driver "INTUN" sets up the necessary input parameters for the subprogram Intun and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Intun", disc 1 

• Calling Syntax 

CALL Intun (N, Eps, X(*), Y(*), Int) 

• Input Parameters 

N Number of data points; N must be at least 3. 

Eps Epsilon tolerance for the solution of the simultaneous equa- 

tions used in the calculation of the integral. 

X(*) Vector containing the domain values of the data points sub- 

scripted from 1 to N; domain values must be in increasing 
order. 

Y(*) Vector containing the range values of the data points sub- 

scripted from 1 to N. 

• Output Parameters 

Int Value of the integral. 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Intun." 
N = Esp = 

X_(l)= X_(N) = 

The data may be corrected from the keyboard (e.g., N = 21, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• For a derivation of the equations and the flow chart upon which this program is based, 
see Greville's article in reference 1. Reference 2 explains the theory of spline functions in 
general with applications in various other areas. 

• Vectors X(*) and Y(*) must be dimensioned in the calling program: 
DIM X(1:N), Y(1:N) where N is the number of data points to be used. 

• The number of data points must be greater than 2, otherwise an error will occur in the 
dimension statement in line XXX DIM B(2:N - 1), S(N), G(2:N - 1) of the subprogram. 

Methods and Formulae 

This subprogram approximates |"f(x)dx where f(x) is represented by discrete functional values 
for unequally-spaced domain values x over the interval [a,b]. 

The method implemented involves fitting a curve through the data points and integrating that 
curve. The curve used is the cubic natural spline function which derives its name from a 
draftsman's mechanical spline. If the spline is considered as a function represented by s(x), 
the second derivative s"(x) approximates the curvature. For the curve through data points 
(xi,yi), (x 2 ,y 2 ), ■••, (Xn,Vn) we want /^ (s"(x)) 2 dx to be minimized in order to achieve the 
"smoothest" curve. 

The spline function with minimum curvature has cubic polynomials between adjacent data 
points. Adjacent polynomials are joined continuously with continuous first and second deriva- 
tives as well as s"(x x ) = s"(x n ) = 0. 

The procedure to determine s(x) involves the iterative solution of a set of simultaneous linear 
equations by the method of successive overrelaxation. You can specify accuracy to which 
these equations are solved. For a detailed discussion of the algorithm, see reference 2. 

jr ,+1 s(x)dx«"£ 1 -|-(x I + 1 -x 1 ){y 1 + y 1 + 1 )- ± (x i + 1 ) 3 s"( Xi )+ s"(x i + 1 ) 
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References 

1. Ralston and Wilf, Mathematical Methods for Digital Computers, Vol. II, New York: John 
Wiley and Sons, 1967, pp. 156-158. 

2. Greville, T.N.E., Ed., "Proceedings of An Advanced Seminar Conducted by the 
Mathematics Research Center", U.S. Army, University of Wisconsin, Madison, October 
7-9, 1968, Theory and Application of Spline Functions, New York, London: Academic 
Press, 1969, pp. 156-167. " 

User Instructions 

1. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "INTUN" and press RUN. 

2. You will be asked to supply entries to the following items: 

- number of data points to be used in the computation 

- error tolerance desired in solving the system of equations 

- x (in increasing order) and y coordinate of each data point. 
Press CONTINUE after each entry. 

3. The program will then print the endpoints of the interval of integration and the value of 
the integral. 
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Examples 

User entries: 



Result: 



User entries: 



Result: 



* OF DATA POINTS-- 5 
ERROR TOLERANCE* i . E-6 



DATA POINTS: [DOMAIN VALUES MUST BE 
IN INCREASING ORDER] 

X( i)~ Q.00000QE+00 Y( i)= . 00 OE-M) 

X ( 2 ) ;::: i. . E + Y ( 2 ) ~ 1.0000 E + 

X( 3)« 2.000000E+00 Y( 3>= 2.000000E+00 

X ( 4 ) ■■■■■■■ 3.00 E + 00 Y ( 4 ) ~ 3 .000 E + 

X ( S ) ::= 4.00000 E + 00 Y ( 5 ) = 4 .0000 E + 

I N T E ORAL F ROM 0.0 T O 4.0 I S 8 . E + 



OF DATA POINTS^ 8 
RROR TOLERANCE-- i . 



DATA POINTS: [DOMAIN VALUES MUST B 
IN INCREASING ORDER 3 



A '•. 



1 ) :••:: 0.00 E + Y (. i ) : " . E + 

X( 2)- 2 . E ■-■ i Y( 2>™ 4 . 0O0OOOE-02 

X( 3):::: 6 . E ■•• 1 Y( 3 ) :::: 3 60 OOOOE-Oi. 

X< 4): ::: .i OOOOOOEMiO Y< A)-" i . 0E-* 

X ( S ) "" i . :i. E •*■ Y < S ) -" 1.21 E •*• 

X< k )'■■■ 1 . SOOOOOE-t 00 Y( h )■■" 2. . 23 0E-M) 

X< 7)- 1 . 6 E ♦■ Y \ 7 ) ~~ 2. S60000E-M1Q 

X ( 8 ) =" : 2 . E •♦• Y ( S ) - 4.00 E ■'<■ 



I N T E C R A i... F R O M 0.0 1 O 2.00 1 8 2.66997 6 E. + 
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Chapter o 
Ordinary Differential Equations 



Introduction 

Description 

A differential equation is a relation among an independent variable and a dependent variable 
and its derivatives. For example, y 1 + xy 2 = 0. The solution is a function y of x satisfying this 
equation and passing through a given point, (x ,Vo)- 

This section contains two methods for solving systems of first order differential equations. 
(Higher order equations may be solved by reducing to a system of first order equations.) 

Routines 

• Kutta - a good, stable method when accuracy requirements are not high or when the 
equation does not involve a large number of calculations. 

• Adams - a method which is more accurate and requires fewer computations but may not 
be as stable in some cases. 
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Runge-Kutta Method 

Description 

This subprogram is a locally fifth-order Runge-Kutta procedure for solving systems of first- 
order ordinary differential equations. 

You must supply the dimension of the system of differential equations, the beginning and end 
points of integration, the initial value vector and the maximum number of integration steps. 



Note 

Before running this subprogram, you must also supply the subpro- 
gram Func containing the user-defined system of equations to be 
solved. If the included driver is used, subprogram Func must be 
stored on the default mass storage device in file "Func". 



Program Usage 

Driver Utilization 

• File Name - "KUTTA", disc 1 

The driver "KUTTA" sets up the necessary input parameters for the subprogram Kutta and 
prints out the resulting outputs. 



Note 

The driver prints every tenth value of the computed vectors. If a 
different number of values is required, line 420 may be changed, 
replacing the 10 by an appropriate number. 



Subprogram Utilization 

• File Name - "Kutta", disc 1 

• Calling Syntax 

CALL Kutta (Idm, A, H, B, Maxstp, Ynt(*), Y(*)) 

• Input Parameters 

Idm Dimension of the system of differential equations. 

A Integration starting point. 

H Integration step size. 

B Integration end point. 

Maxstp Maximum number of integration steps. 

Ynt(*) Initial value vector subscripted from 1 to Idm. 
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• Output Parameters 

Y(*) Y(I,N) is the solution of Ith component evaluated at 

X = A + (N - 1)*H; Y(*) needs to be dimensioned DIM 
Y(Idm, Nb) where Nb is the number of points of integration. 

• Subprograms Required 

Func(Ysv(*), X, Idm, F(*)) Contains user-defined system of equations; used to gener- 
ate new function values, F(*), from old function values, 
Ysv(*); see the Special Considerations section for further 
details. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a check for bad data. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Kutta." 

A = B = 

H = Maxstp = 

The data may be corrected from the keyboard (e.g., H = .01, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• Kutta is a good choice when accuracy requirements are not too high and function eval- 
uations are simple. Kutta requires four function evaluations per integration step, com- 
pared with only two using the Adams predictor-corrector method. 

• The step size, H must be chosen with some care. On the one hand, a small H keeps the 
error due to the Runge-Kutta formulae small. On the other hand, the smaller H is taken, 
the more integration steps we shall have to perform, and the greater the round off error 
is likely to be. 

• Since there may be some cumulative round off error in stepping from end points, A to B 
(i.e., X = A + (N - 1)H), B may not be reached exactly. So, the subprogram is exited 
if X > B - Eps where Eps = IE - 6. If X does not satisfy this condition, the subpro- 
gram will exit when the maximum number of steps have been performed. 

• To define the system of equations to be solved: 

1. PAUSE the current program (if one is running). 

2. Type SCRATCH and press EXECUTE. 

3. In EDIT mode, type in the subprogram Func, pressing ENTER after each line. The 
Subprogram Utilization section explains the parameters used and the examples fol- 
lowing explain how to set up a system of equations in the subprogram. 

4. RE-STORE the subprogram on file "Func". 

5. LOAD file "KUTTA" and press RUN. 

(See "Redefining User Functions" in the Introduction, following the Table of Contents, 
for another way to modify the system of equations.) 
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• A higher order differential equation may be replaced by a system of 1st order equations. 
Assume the equation is of the form: 



d m y dy d 2 y d m_1 y 

— = f(x, y, 



dx m dx dx 2 dx m_1 



.,, . . . , , dy d m - : y 

with the given initial conditions y(x ), — (x ), ..., (x ) 

dx dx" 1 " 1 



We may write the system as follows: 

yi' = fi(x, y l5 y 2 — , y m ) 
y 2 ' = f 2 (x, y : , y 2 — , y m ) 
y m = i m (x, yi, y 2 — , y m ) 

Example 1. 

Compute the solution of van der Pol's equation y" -(.1)(1 - y 2 )y' + y = 
with initial conditions y(0) = 1, y'(0) = 0. 
LetZ = y',thenZ = y" = (.1)(1 - y 2 )y' - y. 

In terms of our subprogram Func, y = Ysv(l) and y' = Ysv(2). 

So, we would code Func as follows: 

SUB Func (Ysv(*),X,Idm,F(*)) 

F: F(l) = Ysv(2) 

F(2) = (.1)*(1 - Ysv(l)*Ysv(l))*Ysv(2) - Ysv(l) 

SUBEND 

Example 2. 

Compute the solution of y" + xy 2 = with initial conditions y(0) = 0, y'(0) = 1. 

LetZ = y' thenZ - y" = -xy 2 . 

In terms of our subprogram Func, y = Ysv(l) and y' = Ysv(2). 

So, we would code Func as follows: 
SUB Func (Ysv(*),X,Idm,F(*)) 
F: F(l) = Ysv(2) 
F(2) = X*Ysv(l)*Ysv(l) 
SUBEND 
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Example 3. 

Compute the solution of 

y'" + 17y" - 10y' + y = with initial conditions y(0) = y'(0) = y"(0) = 0. 

LetZ = y'.thenZ' = y". 

LetW = Z' = y",thenW = Z" = y'" = -17y" + 10y' - y. 

In terms of our subprogram Func, y = Ysv(l), y' = Ysv(2) and y" = Ysv(3). 

So, we would code Func as follows: 

SUB Func (Ysv(*),X,Idm,F(*)) 

F: F(l) - Ysv(2) 

F(2) = Ysv(3) 

F(3) - 17*Ysv(3) + 10*Ysv(2) + Ysv(l) 

SUBEND 
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Methods and Formulae 

This subprogram solves the following system of simultaneous first-order ordinary differential 
equations: 

dyj 

— = f : (x, y 2 , y 2 , ..., y n ) 

dx 

dy 2 

— - =f 2 (x, yj, y 2 , ..., y n ) 

dx 



— - = in(x, yi, y 2 , ..., y n ) 
dx 



with you specifying the initial conditions, x , y^Xo), y 2 (x ), ..., y n (x )- 
The program may also be used to solve an equation of the form: 

d m y dy d 2 y d"-^ 
= f(x,y,— 1—1.,..., 1) 

dx m dx dx 2 dx 1 "- 1 

, dy d m_1 y 

with given initial conditions y(x ), (x ), ..., -(x ) 

dx dx" 1 " 1 

by rewriting the equation as a system of first-order equations as above. A method for doing 
this is provided in the Special Considerations and Programming Hints section. 

Runge-Kutta methods attempt to obtain greater accuracy, and at the same time avoid the 
need for higher derivatives, by evaluating the function f at selected points on each subinterval. 

For the equation y' = f(x,y), y(x ) = y and step size h, approximations y n to y(x + nh) for 
n = 0, 1,2, ... are generated using the recursion formula 

V"n + i = y n + g (ki + 2k 2 + 2k 3 + k 4 ) 
where 

ki = hf(x n , y n ) 

k 2 = hf(x n + 2- y n + 2 k i) 

L -j 

k 3 = hf(x n + 2, y n + -k 2 ) 
k 4 = hf(x n + h, y n + k 3 ) 
The local truncation error is 0(h 5 ). 
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This is a single step method. It requires only the value of y at a point x = x n to find y and y' at 

x = x n + l- 

The above formulae may be generalized to any system of 1st order differential equations. 
Further details may be obtained from the references. 

References 

1. Ralston and Wilf, Mathematical Methods for Digital Computers , Vol. 1, New York: John 
Wiley and Sons, Inc., 1960, p. 115. 

2. Carnahan, Luther, Wilkes, Applied Numerical Methods , New York: John Wiley and 
Sons, Inc., 1969, pp. 363-366. 

User Instructions 

1. If you have not already defined the system of equations you wish to solve, you should do 
so at this time and save on file "Func". 

2. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "KUTTA" and press RUN. 

3. You will be asked to supply entries for the following items: 

- dimension of the system of differential equations 

- lower bound of the integral 

- upper bound of the integral 

- step size for integrating 

- maximum number of iteration steps permitted 

- initial values for y for each dimension of the system. 

4. The program will print the domain values x as well as the values of the integrated 
function at x. 
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Examples 

• Solve y' - xy 1/3 = for 1«= x s= 3, step size = 0.01 subject to the initial condition 
y(l) = 1. Print every tenth value. With y = y sv (l), subprogram Func is set up as 

10 SUBFunc(Ysv(*),X,Idm,F(*)) 
20 F: F(1) = X*Ysv(1)a(1/3) 
30 SUBEND 



User entries: 



I) I 



n F 



YE TEH OF" I) 



LOWER BOUND- 1 
UPPER BOUND"- 3 
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MAX # !'" I M T E: Q E A T 1 N G "f IE P S 



Results: 
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Solve y" + 2yy' = for =s x «= 2, step size = .01, subject to the initial conditions 
y(0) = 0, y'(0) = 1. 

The solution may be obtained by first reducing the equation to a set of simultaneous 
equations using the following substitutions: 

LetZ = y'.thenZ' = y" = -2yy'. 

With y = Ysv(l) and y' = Ysv(2), the subprogram Func is set up as follows: 

10 SUB Func(Ysv(*),X,Idm,F(*)) 

20 F: F(l) = Ysv(2) 
30 F(2)=-2*Ysv(l)*Ysv(2) 

40 SUBEND 



User entries: 



DI HENS I ON OF SYSTEM OF D.E. 

LOWER BOUND™ 

UPPER BOUND™ 2 

STEP SIZE* . Oi 

MAX * OF INTEGRATION STEPS™ 



20 i 



INITIAL VALUES 



Y< i )■"■ . 0O00O0E-* 
Y< 2 )■■■■■■■ i . 0O0 0E* 



Results: 



1)0 GO I. ''-0 
u o (.. i.' (i ()("-• n : 
(ii,i OE -Oi 
. 00000 OE-0 i 
i. i i j u i) h - i 

. 0E- i 

. 00000 0E- :i. 
E - i 

. 000 00 OE-0 1 
E - U 1 
!"i [) : i I ( ! '- 
1 UU 0L+0 
7 Q (i Or f()0 

.300 E ■♦ 
40 E+0 

. SO OE-i 

f;j] |j!i|; E ■■>■ 

7 u i ■♦• U 

Hi! f: ■■: (j {!.['■*■ 

. 9000O0EH5O 

(.! I i) (.!(, (Ih +0 



!} (I i! [) !![: + 

96679? E -02 
973753E-01 
9131. 2 6 E - i 
79949 017 -Ox 
62il72E-0i 
37 496E-01 
04367 BE -01 
64 368 E- 01 
162979E-01 
6iS942E-0i 
499 0E-01 
336S46E •01 
617232E-01 
8S3516E -01 
0Si483E-(n. 
216686E -01 
354 091E-01 
463 06 0E- 01 
S6237SE-01 
64 776E -01 
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• Solve yy" + 3(y') 2 = for =£ x ^ .2, step size = .01, subject to the initial conditions 
y(0) = l,y'(0) = V 4 . 

The solution may be obtained by first reducing the equation to a set of simultaneous 
equations using the following substitutions: 

LetZ =y\ thenZ = y" = -3(y') 2 y. 

With y = Ysv(l) and y' = Ysv(2), the subprogram Func is set up as follows: 

10 SUBFunc(Ysv(*),X,Idm,F(*)) 
20 F: F(l) = Ysv(2) 
30 F(2)=-3*Ysv(2)*Ysv(2)/Ysv(l) 

40 SUBEND 



User entries: 



i) :i: hems :i: on of system of d.e 

LOWER BOUND^ 
UPPER BOUND^ 2 



Results: 



MAX t OF INTEGRATION STEPS^ 2 0:1 



INITIAL VALUES: 

Y ■. i "> ■■■■■ i . [■" -t 

Y ( 2 > •■"• 2 . S it E "■• 1 



. ODOOOOE-v Of) 

1 (I o (i i) I: - ('i 1 

2 CO) E •- (J :'i. 
3. (JOE -01 
4 . 00000 0E- 01 
':':':. . 0E- i 
6 . ODOOOOE-Oi 

7 . 00000 0E- (J 1 

8 . OOOOOOE •■() :'i. 

9 . 0000O0E-Oi 

i . ooooooe+oo 

:'i. . 100 0O0E+0O 
i. . 200 00 0EH)0 
i . 30000 ()E+ 
i . 400000 E +0 
i . 5000 0E+ 
i . 600 00 0E+0 
1 7000 0E+0 
i . 800 0E+0 
:'i. . 90 0E+00 
2 . O0 00OOE+0O 



it ii Oi : "f ii 
24114E 



4 - rA 

06779 0E 

877S7E 

i 6682 E 

i 24683E 

14i8B8E 

1B8292E 

1740SSE 

1B9207E 

20380 IE 

2I7883E 

231493E 

24 4 666 E 

257433E 

2698231": 

28i86iE 

293S69E 

304967E 

3160 7 4E 
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+ 
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+ 

+ 

+ 
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• Solve y" - (.1)(1 - y 2 )y' + y = for « x =£ .2, step size = .01, subject to the initial 
conditions y(0) = 1, y'(0) = 0. 

The solution may be obtained by first reducing the equation to a set of simultaneous 
equations using the following substitutions: 

Let Z = y\ then Z' = y" = (.1)(1 - y 2 )y' - y. With y = Ysv(l) and y' = Ysv(2), 
the subprogram Func is set up as follows: 

10 SUB Func(Ysv(*),X,Idm,F(*)) 

20 F: F(l)=Ysv(2) 

30 F(2) = .l*(l-Ysv(l)*Ysv(l))*Ysv(2)-Ysv(l) 

40 SUBEND 



User entries: 

Din I 

LOUIE 
UP PI 
STE 

MAX 



■NSION OF SYSTEM OF D.E.= 2 
:.R BOUND = 

;;:r bouni>= 2 
nzF.= .01 

OF INTEGRATION STEPS--- 201 



Results: 



INITIAL VALUES: 

Y( i)"- 1 . 00 ODE-MO) 
Y( 2'i^ . 0000 E •*•()() 



1 , (J . *■ (3 

i. . E ■ 1 

,-' not I; 1 '. -0 i 

3 ooooooE--oi 

4 no u fi ,..1 ■' -0 1. 

'■ : (! 0!" •• 1. 

{-> . (i (: E - 1 

.,' (i fi I- (s U: -0 :> 

8 n (i e - 1 

- u f - X 

■\ i ii ii !'i )'. Mi 

1 1 i (Mi 01':. -f 

1 .'■ :> ' . .: I !■- 
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1 / I: ''■ 

:i . GO (O)0 OEM) 

i Mi ;.; n u (jf { u 

= 0(00.0)0 0+00 



g 
a 



ft i\ -3 i.\ '■■ 
950 4 IE™ Gi 
B0 0650L - 01 
55324 6E 01 
21010.90 - 01 
77436 EMU 
2498 09E-01 
641 3E-0 1 
9 3313 7E -01 
1.9 '0)40-1:: -0 1 
364177E -01 
47660 0E- 01 
336993E -01 
S03t:MlE-01 
335432E -01 
918325E -02 
6714900-0 9 

631 023E -01 
68891-02 -01 
M.^OOH. -01 
741940E 01 
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Adams-Bashford-Moulton Method 

Description 

This subprogram is a locally fifth-order Adams-Bashford procedure for solving systems of 
first-order ordinary differential equations. Kutta, a Runge-Kutta procedure, is used as a start- 
er. An optional Adams-Moulton corrector is also included. 

You are required to supply the dimension of the system of differential equations, the begin- 
ning and end points of the integration, the initial value vector and the maximum number of 
integration steps. 



Note 

Before running this subprogram, you must also supply the subpro- 
gram Func containing the user-defined system of equations to be 
solved. If the included driver is used, subprogram Func must be 
stored on the default mass storage device in file "Func". 



Program Usage 

Driver Utilization 

• File Name - "ADAMS", disc 1 

The driver "ADAMS" sets up the necessary input parameters for the subprogram Adams and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Adams", disc 1 

• Calling Syntax 

CALL Adams (Idm, A, H, B, Maxstp, Crktr, Ynt(*), Y(*), Er(*)) 

• Input Parameters 

Idm Dimension of the system of differential equations. 

A Integration starting point. 

H Integration step size. 

B Integration end point. 

Maxstp Maximum number of integration steps. 

Crktr Corrector flag: 

Crktr = 1 implies use corrector. 

Crktr = - 1 implies do not use corrector. 

Ynt(*) Initial value vector, subscripted from 1 to Idm. 

• Output Parameters 

Y(*) Y(I,N) is the solution of Ith component evaluated at 

X = A + (N - 1)*H; Y(*) needs to be dimensioned DIM 
Y(Idm, Nb) where Nb is the number of points of integration. 

Er(*) Error estimate for Y(I,N); Er(*) is subscripted from 1 to Nb 

where Nb is the number of points. 



Ordinary Differential Equations 3-13 



• Subprograms Required 

Kutta (Idm, A, H, B, 

Maxstp (Ynt(*), Y(*)) Used as a starter for subprogram Adams. 

Func (Ysv(*), X, Idm, F(*)) Contains user-defined system of equations; used to gener- 
ate new function values F(*), from old function values, 
Ysv(*). 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Adams." 

A = B = 

H = Crktr = Maxstp = 

The data may be corrected from the keyboard (e.g., H = .01, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• When functional evaluations are costly, Adams is a good choice since only two evalua- 
tions are made per integration step. 

• The step size, H, must be chosen with some care. On the one hand, a small H keeps the 
error due to the Adams-Bashford formulae small. On the other hand, the smaller H is 
taken, the more integration steps we shall have to perform, and the greater the rounding 
error is likely to be. 

• Since there may be some cumulative round off error in stepping from end points A to B 
(i.e.. X = A + (N - 1)*H), B may not be reached exactly. So, the subprogram is exited 
ifX>B-lE-8. 

• The Moulton corrector term is an available option. In most of the test cases tried, one or 
two digits of accuracy were gained. The cost is an additional 50 - 75% in running time. 

• To define the system of equations to be solved: 

1. PAUSE the current program (if one is running). 

2. Type SCRATCH and press EXECUTE. 

3. In EDIT mode, type in the subprogram Func, pressing ENTER after each line. The 
Subprogram Utilization section explains the parameters used and the examples on 
pages 3-8 through 3-11 explain how to set up a system of equations in the sub- 
program. 

4. RE-STORE the subprogram on file "Func". 

5. LOAD file "ADAMS" and press RUN. 

(See "Redefining User Functions" in the Introduction, following the Table of Contents, for 
another way to modify the system of equations. ) 
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Methods and Formulae 

This subprogram solves the following system of simultaneous first-order ordinary differential 
equations: 



^L = f 1 (x,y 1 ,y 2 ,...,y n ) 
dx 

^2_ = f2(x >yil y 2 ,...,y n ) 

dx 



^n_ = f„(x.y 1 ,y 2 ,... I y n ) 
dx 

with you specifying the initial conditions, x , y^xo), y 2 (x ), • •■, y n (x ). 

The program may also be used to solve an equation of the form: 

d m y dy d 2 y d m "V 
= f(x, y, , ,..., ) 

dx m dx dx 2 dx ml 

with given initial conditions 

dy d m_1 y 

y(x ), (xo), ..., -(x ) 

dx dx m_1 



by rewriting the equation as a system of first-order equations as above. A method for doing 
this is provided in the Special Considerations and Programming Hints section. 

For the equation y' = f(x,y), y(x ) = y and step size, h, approximations y n to y(x + nh) for 
n = 0, 1, 2, ... are generated using the recursion formula 

Vn + i = y n + ^ (55f n - 59f n - 1 + 37f n - 2 - 9f n - 3) 

where f, = f(x h y,). 

The local truncation error is 0(h 5 ). 

This is a multistep method; thus, it is not self-starting. We must have four successive values of 
f(x, y) at equally-spaced points before this formula can be used. The subprogram Kutta, a 
Runge-Kutta procedure, is used in Adams to generate these points. 

The above formulae may be generalized to any system of first-order differential equations. 
Further details may be obtained from the references given on page 3-7. 
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User Instructions 

1. If you have not already defined the function you wish to solve, you should do so at this 
time. 

2. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "ADAMS" and press RUN. 

3. You will be asked to supply entries for the following items: 

- dimension of the system of differential equations 

- lower bound of the interval under consideration 

- upper bound of the interval under consideration 

- step size or distance between each point of integration 

- maximum number of integration steps permitted 

- initial values for y for each dimension of the system. 
You may also choose to use the Moulton corrector. 
Press CONTINUE after each entry. 

4. The program will print the domain values x, as well as the values of the integrated 
function at x. 
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Examples 

• See example on p. 3-8 for a statement of the problem. Subprogram Func is set up as 
follows: 

10 SUB Func(Ysv(*),X,Idm,F(*)) 
20 F: F(1) = X*Ysv(1)a(1/3) 
30 SUBEND 

User entries: 

DIMENSION OF SYSTEM OF I) . E . - 1 

LOWER B0UND= 1 

UPPER B0UND= 3 

STEP SIZE> . Oi 

MAX * OF INTEGRATION STEPS™ 201 



INITIAL VALUES: 

Y< i)= i.OO00O0E+00 



Results: 



X Y 

i. OOOOOOE+00 1. OOOOOOE + 00 

1 .iOOOOOE+00 i . i068!7E+00 
1.200000E+00 i.227880E+00 
1 . 300000E+00 i . 364136E+00 
i .400000 E+ 00 i.5i656SE+0 
i . 50000 0E+0 i . 68617.1.E + 
1.600000 E + 00 i . 873982E + 
i . 700000E+00 2. 031045E+00 
i.BOOOOOE + 00 2.308421E+00 
i 900000 E+ 2 . 557 i 87E + 

2 . E+ 2 . 8 2 8427 E + 
2. i00OOOE+O0 3. 123239E+00 
2.200000E+00 3. 442725E+00 
2 .3000 0E+0 3 . 787995E + 
2 . 400000E+00 4 . 160166E+00 
2 . 5 E+ 4.560 359E+ 
2. 600000E+00 4. 989698E+00 
2.700000 E + 5 . 4493 i 2E + 
2 . 800000E+00 5 . 940333E+00 
2 . 90 0E+0 6 . 463894E+0 
3. OOOOOOE+00 7.021132E+00 
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• See example on p. 3-9 for a statement of the problem. Subprogram Func is set up as 
follows: 

10 SUB Func(Ysv(*),X,Idm,-F(*)) 
20 F: F(l)=Ysv(2) 
30 F(2)=-2*Ysv(l)*Ysv(2) 

40 SUBEND 



User entries. DIMENSI0N 0F SY STEM OF D . E . - 2 

LOWER BOUNDS 
UPPER BOUND-::: 2 
STEP SIZE*: .01 
MAX * OF INTEGRATION STEPS* 20i 



INITIAL VALUES-. 

Y( i)= 0.0 00 OE+0 
Y( 2)= i.OOOOOOE+00 

Results: 



. 


000 000E+0 





.00 000 OE+0 


i 


.QOOOOOE-Oi 


9. 


966799E- 


-02 


2. 


OOOOOOE-Oi 


i. 


9737S3E- 


- i 


3 


. OOOOOOE-Oi 


2. 


913126E- 


-01 


A. 


OOOOOOE-Oi 


3. 


799490E- 


•Oi 


5 


. OOOOOOE-Oi 


4. 


.62ii72E- 


-01 


6. 


OOOOOOE-Oi 


5. 


370 496E- 


-Oi 


7 


.000 0OOE-Oi 


6 


. 43678E- 


-01 


8. 


OOOOOOE-Oi 


6. 


640368E- 


-Oi 


9 


.000OO0E-Oi 


7 


. 162979E- 


-01 


i . 


. 00O00OE+0O 


7. 


61S942E- 


■ i 


i 


. 100000E+00 


8 


. 04991E- 


-01 


i. 


. 200000E+00 


8. 


336S47E- 


- i 


i 


. 3000 OE+0 


8 


. 617232E- 


-Oi 


i. 


. 40 000 0E+0 


8. 


. 853S17E- 


■01 


i 


.50000 OE+0 


9 


, 051483E- 


•01 


1 


.600 OE+0 


9. 


2166S6E 


-01 


i 


. 700000E+00 


9 


.354091E- 


-Oi 


i 


.80 000 OE+0 


9. 


468 061 E- 


- .1. 


i 


.900000E+00 


9 


.562375E- 


-Oi 


2 


, 00 000 0E+0 


9. 


.640276E- 


•01 
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• See example on p. 3-10 for a statement of the problem. Subprogram Func is set up as 
follows: 

10 SUBFunc(Ysv(*),X,Idm,F(*)) 

20 F: F(l)=Ysv(2) 
30 F(2)=-3*Ysv(2)*Ysv(2)/Ysv(l) 

40 SUBEND 

User entries: 

DIMENSION OF SYSTEM OF D.E.= 2 

LOWER BOUND- 

UPPER BOUND- 2 

STEP SIZE- . 01 

MAX t OF INTEGRATION STEPS- 20 i 



INITIAL VALUES: 

Y( i>- 1.0 000 00E+00 
Y( 2)= 2.500000E-01 

Results: 






. OOOOOOE + 


1 


.OOOOOOE+0 


i 


. 00000 OE-Oi 


i 


. 024U4E + 00 


2 


. 00 00 OE-Oi 


i 


. 04663SE + 


3 


. 00 000 OE-Oi 


i 


. 067790E+0 


4 


. 0000O0E-O1 


i 


. 0877S7E+0 


5 


. OOOOOOE-Oi 


i 


. 106682E+0 


6 


. OOOOOOE-Oi 


i 


. .i.24683i;::+ oo 


7 


. OOOOOOE-Oi 


1 


. 1418S8E+0 


8 


. OOOOOOE-Oi 


i 


. 158292E+00 


9 


.000000E-0i 


i 


. 1740SSE+0 


i 


. 000000E+00 


i. 


189207E+0 


i 


. iOOQQOE+0 


i 


.203801 E +00 


i . 


200000E+00 


i. 


217883E+0 


i 


. 300 00 0E+0 


i 


23.1493E+0 


:i. . 


400000E+00 


i. 


244666E+0 


1 . 


. S000 0E+0 


1 . 


257433E+0 


i . 


600 00 0E+0 


i. 


269823E+0 


i . 


700 0E+0 


i . 


281861E+00 


i . 


800 00 0E+0 


i . 


293569E+0 


i . 


900000E+00 


i . 


304967E+0 


2. 


O00 000E+OO 


i . 


316 074E+0 
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• See example on p. 3-11 for a statement of the problem. Subprogram Func is set up as 
follows: 

10 SUB Func(Ysv(*),X,Idm,F(*)) 

20 F: F(l)=Ysv(2) 

30 F(2) = .1*(1 - Ysv(l)*Ysv(l))*Ysv(2) - Ysv(l) 

40 SUBEND 



User entries: 



DIMENSION OF SYSTEM OF D.E.= 2 

LOWER BOUND" 

UPPER BOUND* 2 

STEP SIZE= .01 

MAX * OF INTEGRATION STEPS" 20 i 



INITIAL VALUES 



Y< i>» 

Y( 2) = 



1.0 00 0E+0 
O.OOOOOOE+00 



Results: 




1 

2 
3 
4 
5 
6 
7 



00 

00 



00 
00 




8. 
00 




9 
i 
1. 
i 
i 
1 
1 
i 
j. 
i 
i 



10 
2 
30 
4 
SO 
60 
70 
80 
90 




0O00E+00 
0000E-01 
0000E-01 
000 0E-01 
0O0OE-O1 
00 0OE-01 
00OE-01 
OOOOE-Oi 
00 0E-01 
00 0E-01 
00 0E+O0 
OOOOE+0 
0OE+O0 
00 0E+0 
OO0OE+0 
OOOOE+0 
OOOOE+0 
00 0E+0 
OOOOE+00 
OOOOE+0 
O000E+OO 



OOOOOOE 
9S0041E 
800650E 
S53246E 
210U9E 
774360E 
249809E 
6410 03E 
9S3137E 
19204SE 
5. 36 4177 E 
4.476600E 
S36993E 
55364 IE 
S35432E 
91832SE 
671 498 E 
63102SE 
68891 2E 
72965 0E 



741948E- 



+ 00 
-01 
-01 
-01 
•••■ 1 
-01 

•01 
-01 
•••■ 1 
-01 
-01 
-01 

•• .1. 
-01 
-01 
-02 
-02 
-01 
-01 
-01 
01 
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Chapter 4 
Linear Algebraic Systems 



Introduction 

Description 

This section contains a wide variety of routines dealing with matrices and systems of simul- 
taneous linear equations. Significant savings in storage are obtained in the matrix inversion 
routines by taking advantage of special types of matrices, such as positive definite. An addi- 
tional digit or two of accuracy may also be obtained in many cases. There are also routines for 
storing symmetric matrices using minimal memory and for transposing a matrix in place, that 
is, replacing a matrix by its transpose. 

Routines 

• Ludsht - solves the system of equations Ax = b using triangular decomposition. 

• Decomp - computes the triangular decomposition of nonsingular matrix A; contained in 
file "Ludsht.". 

• Solve - finds an approximate solution to a single system of equations, Ax = b; contained 
in file "Ludsht". 

• Posdef - given a symmetric, positive definite matrix, A, stored in symmetric storage 
mode, will solve the system of equations Ax = b using Cholesky's Method. 

• Choles - performs Cholesky decomposition on a symmetric, positive definite matrix 
stored in symmetric storage mode; contained in file "Posdef". 

• Solcho - solves a Cholesky system in symmetric storage mode; contained in file 
"Posdef". 

• Pinver - finds the inverse of a positive definite matrix stored in symmetric storage mode 
in vector S(*). 

• Triang - given a lower triangular matrix stored in symmetric storage mode vector S(*), 
this routine will multiply S T S; contained in file "Pinver". 

• Invers - finds the inverse of a lower triangular matrix stored in symmetric storage mode 
vector S(*); contained in file "Pinver". 

• Storag - will store a symmetric matrix, A, in a vector, S, in symmetric storage mode, or a 
vector, S, into a symmetric matrix, A. 
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Linear Equation Solver-Minimum Storage 

Description 

This subprogram solves the system of equations AX = B using triangular decomposition. The 
triangular matrices L and U overwrite A and the solution matrix X overwrites B. This saves a 
significant amount of memory, but in the process, the values contained in matrices A and B 
are destroyed. 

Program Usage 

Driver Utilization 

• File Name - "LUDSHT", disc 1 

The driver "LUDSHT" sets up the necessary input parameters for the subprogram Ludsht 
and prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Ludsht", disc 1 

• Calling Syntax 

Call Ludsht (A(*),B(*),N,M) 

• Input Parameters 



A(*)' 

B(*)' 

N 
M 



Array containing the nonsingular matrix A; dimensioned 
A(1:N, 1:N). 

Array containing the coefficient matrix B; dimensioned 
B(1:N, 1:M). 

Number of rows in B(*). 

Number of columns in B(*). 



• Output Parameters 

M*Y 
B(*Y 



Array containing the Lu decomposition of A(*). 

Array containing the solution matrix X to the system 
AX = B. 



Subprograms Required 
Decomp 

Solve 



Computes the triangular decomposition of A(*) using Gaus- 
sian Elimination. 

Given the triangular decomposition of A(*), Solve finds an 
approximate solution to a single system of equations 
AX = B. 



Further explanation of these subprograms is given on pages 4-6 through 4-10. 



1 A(*) and B(*) are both input and output parameters. 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Ludsht." 
M = N = 

The data may be corrected from the keyboard (e.g., N = 5, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• One way to solve the system of linear equations AX = B is to compute the inverse of A 
and then multiply A" 1 by B. This method may appear especially attractive if several right- 
hand sides B k are involved since the inverse need be computed only once. However, a 
set of linear equation solving procedures - such as Decomp and Solve - can accomplish 
this task with fewer operations and with greater accuracy. Once L and U have been 
computed, the solution of LUX = B requires n 2 - n multiplications and n divisions or a 
total of n 2 multiplicative operations. Moreover, once A" 1 has been computed, the evalua- 
tion of A _1 B requires n 2 multiplicative operations also. Thus both methods require appro- 
ximately the same number of operations at this point. But the initial calculations of 
LU and A" 1 require about in 3 and n 3 multiplicative operations, respectively. 

• This subprogram is designed to save as much storage as possible. But care must be taken 
since both matrix A and matrix B are destroyed in the subprogram. 

Methods and Formulae 

Subprogram Ludsht begins by finding the triangular decomposition of A using subprogram 
Decomp. The decomposition Lu(*) is stored in A(*). 

AX = B is then solved one column at a time using subprogram Solve. The solution vector is 
then restored in B(*), 

Reference 

1. Forsythe, G. and Moler, C, Computer Solution of Linear Algebraic Systems, Engle- 
wood Cliffs, N.J.: Prentice Hall, Inc., Ch. 9, 11. 

User Instructions 

1. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "LUDSHT" and press RUN. 

2. You will be asked to supply entries for the following items: 

- dimensions of the coefficient matrix B (row and column) 

- elements of nonsingular matrix A 

- elements of coefficient matrix B 
Press CONTINUE after each entry. 

3. The program will then print the solution matrix X. 
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Examples 

User entries: 



D [Mi M:: IT) ;•■•■; 0' B< * '■'-" 



.'.'! Y" o 



flATRIX A 



1. 00 9 00!: * -.■'■' i' m ■; <!'.. ■;■ i.'i ':' fj fl Hi K) (i 

:'i. . ', U H F ; - 1) - /' . M.; On I -Ml i Mil Mi 

2 . E i- - ! £ + ,S E - 

1 : ! ; ; i.i (] I v- ;i i MOOM'uMMO '■■ >; jinjii r ■ 'M 

1. MM .- MMH ■ iMl •-•'! 00 00 0;' MM , Oil ooooi <mo 
.3 o i.i fl o n 0E J > n o 



i!- ! 



v . UOI.i !■:. + !} I '' OU-JOlH.il: *'.; : ! ■• -v JO Mi* MM 

I"! . i'l n l; i| J ()! '-Oil ■ . i! if (i !..'' "' !: Ml ' 0!; '• ) 

B . 00 OOE M) . 00 OIEHH) 



Results: 

;i A i !■•! J. 



i. (i E •* ■•- 2 . 8 8 4 6 :i. 3 V-.. ■• ■■ i 1.0 0. •*■ 
'."'' M '.- . : .''i - ;'i. 00 00 0! M10 -• i 'M.-MOM- -■ . ; 

1 . i.!F-* n •■■■;i. i <' ,M'; ■-},-., f. -Oi. 
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User entries: 

DIMENSIONS OF B(*)= 3 X i 



MATRIX A 



Results: 



3.300000 E ■«■ i. i . 6 E + :i. 7.2000 (1 E + Q 1 
•2.400 £ i- :i. ■- i . E + :i. ■■■■ 5.700000 E + :i. 
• 8 . E ■♦■ -• 4.00 E + ■■- i . 7 E + i 



MATRIX B- 

-• 3 . S 9 E + 02 2 . 8 1 E + 2 (3.30 E '<■ :i. 

MATRIX X : 

i . E + •- 2.0000 E + ■••• 3.00 E ■♦■ 
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Subprogram Decorap 

Given a nonsingular matrix A, subprogram Decomp computes the triangular decomposition of 
A using Gaussian Elimination. The triangular matrices L and U are calculated as well as the 
permutation matrix P such that LU = PA. 

Subprogram Utilization 

• FileName - contained in file "Ludsht", disc 1 

• Calling Syntax 

CALL Decomp (N, A(*), Lu(*), Ips(*)) 

• Input Parameters 

N Size of matrix A. 

A(*) Array containing nonsingular matrix A; dimensioned A(1:N, 

1:N). 

• Output Parameters 

Lu(*) Array storing (L-I) and U, the triangular matrices in the 

triangular decomposition; dimensioned Lu(l:N, 1:N). 

Ips(*) Vector containing the permuted indices; subscripted from 1 

to N. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Decomp." 

N = 

The data may be corrected from the keyboard (e.g., N = 10, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• If the subprogram detects a row of zeros, the following error message will be printed and 
the program will pause: 

"ERROR IN SUBPROGRAM Decomp." 
"MATRIX WITH ZERO ROW." 

This test is made at the beginning of the subprogram before any changes are made in 
any of the matrices. It indicates some problem in the calling program. You may want to 
restart the program to correct this. 
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• In the Gaussian Elimination section, there are two tests performed to check that A(*) is 
not a machine singular matrix. If either test fails, the following error message is printed 
and the program pauses: 

"ERROR IN SUBPROGRAM Decomp." 
"MATRIX IS MACHINE SINGULAR." 

The difficulty here is more complex. Certain nonsingular matrices may be made singular 
as a result of the perturbations introduced by round-off error. If this is the case, a more 
specific algorithm may be required to solve the system of equations. See the references 
to Golub and to Golub and Kahan for further details. 

More likely, a truly singular input matrix will be perturbed into a neighboring nonsingular 
matrix by the round-off since normally round-off modifies some element of the pivotal 
column to a non-zero value. 

Methods and Formulae 

Given a nonsingular matrix A, subprogram Decomp computes the triangular decomposition 
using Gaussian Elimination. The triangular matrices L and U are calculated as well as the 
permutation matrix P such that LU = PA. 

Temporarily ignoring scaling and pivoting, the central calculation, the elimination, can be 
expressed by 

FORJ = K+ltoN 

A(I,J) = A(I,J) - (A(I,K)/A(K,K))*A(K,J) 
NEXT J 

This operation is carried out by the innermost FOR-NEXT statement. The multipliers, (A(I,K)/ 
A(K,K)), are saved in the lower triangular matrix L. 

In Decomp, the element of largest absolute value in each row of the matrix is found and its 
reciprocal is recorded in vector Scales(*). But no actual scaling is carried out. Instead, these 
scale factors are used for choosing the pivot element only. This technique has two favorable 
consequences: exact powers of the machine base are not needed for scaling, and the scale 
factors do not have to be applied to the right-hand sides. 

The same type of consideration is involved in pivoting. The array Ips(*) is initialized so that 
Ips(I) = I. 

During the elimination, the largest element in the column is chosen as the pivot element, but 
the rows are not actually interchanged. The corresponding elements of Ips(*) are inter- 
changed instead. We then refer to A(Ips(I),J) instead of A(I,J). This involves no great loss of 
time as long as all inner loops are on the column subscript J. We gain the time that would be 
required to carry out the interchange. 

Finally, L - 1 and U are stored in matrix Lu. 
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References 



1. Golub, G., "Numerical Methods for Solving Linear Least Squares Problems" Numer 
Math ., Vol. 7 (1965) pp. 206-216. 

2. Golub, G. and Kahan, W., "Calculating the Singular Values and Pseudo-inverse of a 
Matrix", J. SIAM Numerical Analysis Series B, Vol. 2 (1965) pp. 205-224. 
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Subprogram Solve 

Given the triangular decomposition of a nonsingular matrix A, stored in matrix Lu, this sub- 
program will find an approximate solution to a single system of equations, AX = B. 

Subprogram Utilization 

• File Name - contained in file "Ludsht", disc 1 

• Calling Syntax 

CALL Solve (N, Lu(*), B(*), X(*), Ips(*)) 

• Input Parameters 

N Order of matrix Lu(*). 

Lu(*) Array containing the triangular decomposition of the non- 

singular matrix A; dimensioned Lu(l:N, 1:N). 

B(*) Vector containing the coefficients B of AX = B; subscripted 

from 1 to N. 
Ips(*) Vector containing the permuted indices from subprogram 

Decomp; subscripted from 1 to N. 

• Output Parameters 

X(*) Vector containing the solution to AX = B; subscripted from 

1 to N. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Solve." 

N = 

The data may be corrected from the keyboard (e.g., N = 10, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• If the subprogram detects a division by zero, the following error message is printed and 
the program will pause: 

"ERROR IN SUBPROGRAM Solve." 
"DIVISION BY ZERO DETECTED." 

This indicates that the matrix is machine singular. Certain nonsingular matrices may be 
made singular as a result of the perturbation introduced by round-off error. If this is the 
case, a more specific algorithm may be required to solve the system of equations. See 
the references to Golub and to Golub and Kahan (page 4-8) for further details. 
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Methods and Formulae 

Solve uses the Lu factorization from Decomp to find an approximate solution to a single 
system of equations, AX = B. 

Solve consists of two steps. The first solves the .lower triangular system LY = B. The second is 
the back solution, i.e., the solution of the upper triangular system UX = Y. The intermediate 
vector Y is stored in X, and the right-hand side B is not altered. 



Linear Algebraic Systems 4-11 



Positive Definite Matrices 

Description 

Given a symmetric, positive definite matrix, A, stored in symmetric storage mode, Posdef will 
solve the system of equations AX = B using Cholesky's Method. The value of A is over- 
written. 

Program Usage 

Driver Utilization 

• File Name - "POSDEF", disc 1 

The driver "POSDEF" sets up the necessary input parameters for the subprogram Posdef and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Posdef", disc 1 

• Calling Syntax 

CALL Posdef (S(*), B(*), R, C) 

• Input Parameters 



Si*) 1 

B(*)' 

R 
C 

< Output Parameters 
S(*Y 

> Subprograms Required 

Choles 
Solcho 



Vector containing symmetric, positive-definite matrix, stored 

in symmetric storage mode; subscripted from 1 to 

(R + l)R/2. 

Array containing coefficient matrix B; dimensioned B(1:R, 

1:C). 

Number of rows of B(*). 

Number of columns of B(*). 



Array containing Cholesky decomposition, S = GG T where 
G is a lower triangular matrix; subscripted from 1 to 
(R + l)R/2. 

Array containing solution matrix X to AX = B; dimensioned 
B(1:R, 1:C). 



Performs Cholesky decomposition on S(*). 

Solves Cholesky system in symmetric storage mode. 



Further explanation of these subprograms is given on pages 4-15 through 4-17. 



1 S(*) and B(*) are both input and output parameters. 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Posdef." 
R = C = 

The data may be corrected from the keyboard (e.g., R = 5, EXECUTE.) When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• Because of the symmetry of the positive definite matrix, it is necessary to store only 
\ (n + 1), or slightly over half, of its elements, resulting in an important saving of storage 
for large matrices. Since some additional multiplications and additions are required to 
manipulate the vector, there is an increase in execution time. 

• Given a symmetric matrix A, the symmetric storage mode vector is formed as follows: 



A(l, 1) 


A(l,2) 


A(l,3) . 


.. A(l,n) 


A(2, 1) 


A(2, 2) 


A(2,3) . 


.. A(2, n) 


A(3, 1) 


A(3, 2) 


A(3, 3) . 


.. A(3, n) 



A(n, 1) A(n, 2) A(n, 3) ... A(n, n) 

i.e., A(l, 1), A(2, 1), A(2, 2), A(3, 1), A(3, 2), A(3, 3), ..., A(n, 1), A(n, n). Subprogram 
Storag may be used to convert a full storage matrix to symmetric storage mode and vice 
versa. 

Methods and Formulae 

Any positive definite matrix A has a unique decomposition in the form A = GG T , where G is 
a lower triangular matrix with positive diagonal elements. The algorithm is: 

FOR J = 1 to N 

G(J, J) = SQR(A(J, J) - 2 (G(J, K)) 2 ) 
FOR I = J + 1 TO N 

G(I, J) = (A(I, J) - 2 G(I, K)*G(J, K))/G(J, J) 

NEXT I 
NEXT J 

The algorithm is Cholesky's method or the square-root method for factoring a positive definite 
matrix. It is stored in subprogram Choles. 

Given a linear system AX = B, where A is a positive definite matrix, subprogram Posdef first 
calls Choles to factor A into GG T . Then subprogram Solcho is called to forward eliminate and 
back substitute to find the solution matrix X. 
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User Instructions 

1. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "POSDEF" and press RUN. 

2. You will be asked to supply entries for the following items: 

- dimensions of coefficient matrix B (row and column) 

- elements of the symmetric storage vector S 

- elements of coefficient matrix B. 
Press CONTINUE after each entry. 

3. The program will then print the solution matrix X where columns of X are the solutions 
to corresponding columns of B. 



Examples 

/ 4 2 -2\ 
Given the positive definite Matrix A = [ 2 10 5 1 

V-2 5 6/ 

= (S3 6 ) , 
\30 / 



and B = I 33 I , solve AX = B. 



User entries: 

DIMENSIONS OF B(*>= 3 X i 



VECTOR S : 

4.000000 E f 00 2 . E + 1 . E + i 
•■•■ 2.00 E + 00 S . E + 6 I) E + 



MATRIX B 



■6. 00OOOOE+0O 3.300000E+01 3. 000000E>0i 



Results: w, TD ., 

■•• 1 .000 E + 00 2.00 E + 3.00 E + 
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/ 1 1/2 1/3 \ 

• Given the positive definite matrix A = I 1/2 1/3 1/4 I 

Vl/3 1/4 1/5/ 



andB = 



/I U u\ 
= I 1 ), 
V 1 / 



solve AX = B. 



User entries: 



DIMENSIONS OF' B(% 



y EC TOR S 



1 filj i'i Q(i j"' f 



> - U L: • i. 
;■: '■■', o u E o i 



333333E-01 
fi (j (i i I) 1 



I X B 



1 () I;.. + 

f.i o n o n oui.-t- D o 
i) n o o j fi e. <■• o o 



(i (i 
:i. t : (i 
(i 1; 



) u (5 o o n fin or -t n o 
)() OiJiMlOOl ' n 
H) i . 00 00 ()()E -f 00 



Results: 



n ! -i 



c 
60( 
00 ( 



061E+0 

32E-H) 1 
3 E -^ 1 



3 60 03;.'!: 
.1. .930 0171 
1 . 8 1 l : 



01 



02 



3.000030 1: 

•i.eoooisi; 

1 B0 014!i 



+ 1 
■><■ 2. 
■*• 2 
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Subprogram Choles 

This subprogram performs Cholesky decomposition on a symmetric, positive definite matrix 
stored in symmetric storage mode. 

Subprogram Utilization 

• File Name - contained in file "Posdef", disc 1 

• Calling Syntax 
CALL Choles (G(*), R) 

• Input Parameters 

G(*Y Vector containing the symmetric, positive definite matrix 

stored in symmetric storage mode; subscripted from 1 to 
(R + l)R/2. 

R Order of positive definite matrix. 

• Output Parameters 

G(*) ! Vector containing the lower triangular matrix G where GG T 

is the decomposition of the positive definite matrix. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Choles." 
R = 

The data may be corrected from the keyboard (e.g., R = 7, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• If the subprogram detects a nonpositive diagonal element, the following error message 
will be printed and the program pauses: 

"ERROR IN SUBPROGRAM Choles." 

"MATRIX IS NOT MACHINE POSITIVE DEFINITE." 

This is a fatal error. There is no recovery. The PAUSE was inserted to enable you to 
check values of the variables in the subprogram. 

• Given a symmetric matrix A, the symmetric storage mode vector is formed as follows: 

A(l, 1), A(2, 1), A(2, 2), A(3, 1), A(3, 2), A(3, 3), ..., A(n, 1), ..., A(n, n) 



1 G(*) is both an input and output parameter. 
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Methods and Formulae 

Any positive definite matrix A has a unique decomposition in the form A = GG T , where G is 
a lower triangular matrix with positive diagonal elements. The algorithm is: 

FOR J = 1 to N 

j 
G(J, J) = SQR(A(J, J) - 2 (G(J, K)) 2 ) 

K=l 

FOR I = J + 1 TO N 

j-i 
G(I, J) = (A(I, J) - 2 G(I, K)*G(J, K))/G(J, J) 

NEXT I 
NEXT J 

This algorithm is Cholesky's method or the square-root method for factoring a positive defi- 
nite matrix. 
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Subprogram Solcho 

This subprogram solves a Cholesky system in symmetric storage mode; i.e., given a symmet- 
ric, positive definite matrix, A = GG T , stored in symmetric storage mode, solve AX = B. 

Subprogram Utilization 

• File Name - contained in file "Posdef", disc 1 

• Calling Syntax 

CALL Solcho (G(*), B(*), R, C) 

• Input Parameters 

G(*) Vector containing the lower triangular matrix G stored in 

symmetric storage mode, where GG T is the Cholesky de- 
composition of the matrix; subscripted from 1 to (R + 1)R/ 
2. 
B(*) J Array containing the coefficient matrix B; dimensioned 

B(1:R, 1:C). 
R Number of rows in B(*). 

C Number of columns in B(*). 

• Output Parameters 

B(*)' Array containing the solution matrix X of AX = B. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Solcho." 
R = C = 

The data may be corrected from the keyboard (e.g., R = 7, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• The data stored in the coefficient matrix B is destroyed and replaced by the solution 
matrix X. 

• Given a symmetric matrix A, the symmetric storage mode vector is formed as follows: 

A(l, 1), A(2, 1), A(2, 2), A(3, 1), A(3, 2), A(3, 3), ..., A(n, 1), ..., A(n, n). 

Methods and Formulae 

Solcho uses the GG T factorization from subprogram Choles to find an approximate solution to 
the system of equations, AX = B. Both A, stored in symmetric storage mode, and B, the 
coefficient matrix, are destroyed. 

Solcho consists of two steps. The first solves the lower triangular system GY = B. The second 
is the back substitution, i.e., the solution of the upper triangular system G T X = Y. 

1 B(*) is both an input and output parameter. 
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Positive Definite Matrix 

Description 

Pinver finds the inverse of a positive definite matrix stored in symmetric storage mode in 
vector S(*). The inverse overwrites the values contained in S(*). 

This subprogram uses Cholesky's method to decompose the positive definite matrix. 

The advantage of using this algorithm to find the inverse of a positive definite matrix is two- 
fold: 

1. Significant memory is saved. 

2. Since pivoting is not required in the Gaussian elimination, there is less roundoff effect 
on the computed inverse. 

Program Usage 

Driver Utilization 

• File Name - "PINVER", disc 1 

The driver "PINVER" sets up the necessary input parameters for the subprogram Pinver and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Pinver", disc 1 

• Calling Syntax 
CALL Pinver (S(*), R) 

• Input Parameters 



S(*)' 
R 



Vector containing positive definite matrix stored in symmet- 
ric storage mode, subscripted from 1 to (R + l)R/2. 

Order of positive definite matrix; i.e., number of rows or 
number of columns. 



• Output Parameters 
S(*)' 



Vector containing the inverse of the initial positive definite 
matrix stored in symmetric storage mode. 



1 S(t) is both an input and output parameter. 
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• Subprograms Required 

Choles Performs Cholesky decomposition on S(*). 

Invers Finds the inverse of the lower triangular matrix S(*) stored 

in symmetric storage mode. The inverse, also a lower 
triangular matrix, is again stored in S(*). 

Triang Multiplies S T S where S is a lower triangular matrix stored in 

symmetric storage mode. The result is restored in S(*). 

Further explanation of these subprograms is given on pages 4-15 through 4-16 and 4-22 
through 4-23. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Pinver." 

R = 

The data may be corrected from the keyboard (e.g., R = 4, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• This subprogram takes specific advantage of the nature of the matrix. As a result, a 
minimum number of operations are used. Hence, you may expect better accuracy than 
the more general matrix routines. 

• This subprogram also attempts to minimize the storage requirements. The symmetric 
storage vector S(*), containing the positive definite matrix, is repeatedly overwritten. 

• Given a symmetric matrix A, the symmetric storage mode vector is formed as follows: 



A(l, 1) 


A(l,2) 


A(l,3) . 


.. A(l,n) 


A(2, 1) 


A(2, 2) 


A(2, 3) . 


.. A(2, n) 


A(3, 1) 


A(3, 2) 


A(3,3) . 


.. A(3, n) 



A(n, 1) A(n, 2) A(n, 3) ... A(n, n) 

i.e., A(l, 1), A(2, 1), A(2, 2), A(3, 1), A(3, 2), A(3, 3), ..., A(n, 1), ..., A(n, n). Subprog- 
ram Storag (p. 4-24) may be used to convert a full storage matrix to symmetric storage 
mode and vice versa. 



1 S(») is both an input and output parameter. 
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Methods and Formulae 

Any positive definite matrix A, has a unique decomposition in the form A = SS T , where S is a 
lower triangular matrix with positive diagonal elements. This algorithm is Cholesky's method 
or the square-root method for factoring a positive definite matrix. Subprogram Choles is 
based on this algorithm. 

Subprogram Invers is then called to find the inverse of the lower triangular matrix S(*) stored 
in symmetric storage mode. (The symmetric storage mode is used even though the matrix is 
lower triangular and not symmetric.) The inverse, also a lower triangular matrix, is again 
stored in S(*). 

Finally subprogram Triang multiplies S T S to get the inverse of the original matrix: 

A + SS T implies A^ 1 = (SS 7 )" 1 = (S 7 )^ 1 = (S^YS 1 

User Instructions 

1. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "PINVER" and press RUN. 

2. You will be asked to supply entries for the following items: 

- order of positive definite matrix A (number of rows or number of columns) 

- elements of vector S in symmetric storage mode. 
Press CONTINUE after each entry. 

3. The program will print the inverse stored in symmetric storage mode vector S. 
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Examples 



Find the inverse of the following symmetric positive definite matrix: 



/ 1 1/2 1/3 \ 
A = I 1/2 1/3 1/4 I 

\ 1/3 1/4 1/5 / 



User entries: 



ORDER OF A<*> = 3 



VECTOR S: [IN SYMMETRIC STORAGE MODE] 

i.OOOOOOE+00 5.000000E~Oi 3.333333E-01 
3.333333E-01 2 . 500000E-01 2 . O0O00E-O.1. 



Results: 



INVERSE S 



[IN SYMMETRIC STORAGE MODE] 



9. OOOOOOE+00 
3.000000E+01 



■3.600000E+01 
■i.SOOOOOE+02 



1.920000E+02 
i .800 000E+02 



Find the inverse of the following symmetric positive definite matrix: 



( 



1 

1/2 
1/3 
1/4 
1/5 



1/2 
1/3 
1/4 
1/5 
1/6 



1/3 

1/4 
1/5 
1/6 
1/7 



1/4 
1/5 
1/6 
1/7 
1/8 



1/5 
1/6 
1/7 
1/8 
1/9 



) 



User entries: 



ORDER OF A(*> = 5 



VECTOR S 



[IN SYMMETRIC STORAGE MODE] 



1 .OOOOOOE+00 
3.333333E-01 
2.S000O0E-O1 
i . 42857iE-0i 
1 . 428S71E-01 



5. 



OOO0 0OE- 
50 0QE- 
O0OO0OE- 
O0OO0OE- 
250 00E- 



Oi 
01 
01 
0.1. 

01 



3 

2. 

1 

1. 

1 



333333E- 

OOOO00E 

666667E- 

666667E 

iiiiiiE- 



■01 
01 
01 
0:1. 
01 



Results: 



I. IN SYMMETRIC STORAGE MODE] 



2.SD0 000E+01 



i 

1 
1 



050 0E+03 
400000E+03 
792000E+05 



-3.000000E+02 4 

■1.890000E+04 7 

2.688000E+04 -1 

6. 300000 E+ 02 -i 



80 00 0E+03 
93B0 0E+04 
176000E+05 
260000E+04 



5.670000E+04 -8 . 820000E+04 4.410000E+04 
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Subprogram Triang 

Given a lower triangular matrix stored in symmetric storage mode vector S(*), this subpro- 
gram will multiply S T S. The result, a symmetric matrix, will again be restored in vector S(*). 

Subprogram Utilization 

• File Name - contained in file "Pinver", disc 1 

• Calling Syntax 
CALL Triang (S(*), R) 

• Input Parameters 

Sf*) 1 Vector containing lower triangular matrix stored in symmet- 

ric storage mode. 

R Order of full storage mode matrix. 

• Output Parameters 

S(*)' Product of S T S stored in symmetric storage mode. 



Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Triang." 
R - 

The data may be corrected from the keyboard (e.g., R = 5, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

Methods and Formulae 

Given a lower triangular matrix stored in symmetric storage mode vector S(*), Triang multi- 

(z)- 



plies S T S. The result is restored inS( z ). 



FOR I = 1 TO R 
FOR J = I TO R 

R 

S(J*(J - l)/2 + I) = 2. S(L*(L - l)/2 + I)*S(L*(L - l)/2 + J) 

NEXT J 
NEXT I 



L-J 



1 S(*) is both an input and output parameter. 
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Subprogram Invers 

Invers finds the inverse of a lower triangular matrix stored in symmetric storage mode vector 
S(*). The inverse, also a lower triangular matrix, is restored in vector S(*). 

• File Name - contained in file "Pinver", disc 1 

• Calling Syntax 
CALL Invers (S(*),R) 

• Input Parameters 

Si*) 1 Vector containing a lower triangular matrix stored in sym- 

metric storage mode; subscripted from 1 to (R + l)R/2. 

R Order of full storage mode matrix. 

• Output Parameters 

Si*) 1 Vector containing the inverse of the input matrix. The in- 

verse is also a lower triangular matrix and is stored in sym- 
metric storage mode. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Invers." 
R = 

The data may be corrected from the keyboard (e.g., R = 5, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

Methods and Formulae 

Invers finds the inverse of a lower triangular matrix stored in symmetric storage mode vector 
S(*). The inverse, also a lower triangular matrix, is restored in S(*). (The lower triangular 
matrix, of course, is not symmetric. The symmetric storage mode is used solely to save space. ) 

First, the diagonal elements of the inverse are found: 
FOR I = 1 TO R 

S(I*(I - l)/2 + I) = 1/S(I*(I - l)/2 + I) 
NEXT I 

Then, the off-diagonal elements are calculated: 

FOR I = 2 TO R 

FOR J = 1 TO I - 1 

S(I*(I - l)/2 + J) = - S(I*(I - l)/2 + I)* 

i-i 
2 S(I*(I - l)/2 + L)*S(L - l)/2 + J) 

L-J 

NEXT J 
NEXT I 

1 (*) is both an input and output parameter. 
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Symmetric Storage Mode 

Description 

Given a symmetric matrix, A, this subprogram will store A in a vector, S, in symmetric storage 
mode. If A is an n x n matrix, S is a vector with n(n + l)/2 elements. Likewise, given S, this 
subprogram will convert to the symmetric matrix A. 

Program Usage 

Driver Utilization 

• File Name - "STORAG", disc 1 

The driver "STORAG" sets up the necessary input parameters for the subprogram Storag 
and prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Storag", disc 1 

• Calling Syntax 

CALL Storag (A(*), S(*), N, Fig) 

• Input Parameters 

A(*) Array containing symmetric matrix in full storage mode; 

dimensioned A(1:N, 1:N). 

or 

S(*) Array containing vector in symmetric storage mode; dimen- 

sioned S(1:(N + l)N/2). 

N Dimension of symmetric matrix in full storage mode. 

Fig Fig = 1, convert from symmetric to full storage mode. 

Fig = - 1, convert from full to symmetric storage mode. 

• Output Parameters 

A(*) Array containing symmetric matrix in full storage mode. 

or 
S(*) Array containing vector in symmetric storage mode. 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Storag" 

N = Fig = 

The data may be corrected from the keyboard (e.g., N = 6, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• This subprogram saves a significant amount of storage when dealing with large arrays. 
An n x n matrix A requires storage for n 2 elements while vector S contains only 
n(n + l)/2 elements. A large matrix could be stored on a mass storage device in sym- 
metric storage mode and then accessed in particular programs. The cost for the savings 
in storage is time — additional operations are required to convert a particular element in 
the vector S to an element in matrix A. 

Methods and Formulae 

Given the symmetric matrix A and the symmetric storage mode vector S, the A(I, J) element 
of A is stored in S(T) where 

T = I*(I - l)/2 + J if I ss J 
T = J*(J - l)/2 + IifK J 

For example, 



Let A = 1 l 5 6 8 I , then 



(1 2 3 4\ 

2 5 6 8 1 

3 6 7 9 1 

4 8 9 0/ 



S(l) = 


= 1 


S(6) = 7 


S(2) -- 


= 2 


S(7) = 4 


S(3) = 


= 5 


S(8) = 8 


S(4) = 


= 3 


S(9) = 9 


S(5) = 


= 6 


S(10) = 0. 



User Instructions 

1. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "STORAG" and press RUN. 

2. You will be asked to supply entries for the following items: 

- order of the symmetric matrix in full storage mode 

- elements of the symmetric storage mode vector if conversion is from symmetric to 
full storage mode 

- elements of the full storage matrix if the conversion is from full to symmetric storage 
mode. 

Press CONTINUE after each entry. 

3. The program will print the inverted matrix. 
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Examples 

• Full to symmetric storage mode 



User entries: 



ORDER OF MATRIX IN FIR- 



ST OR AGE MOD!:: 



FULL STORAGE MATRIX 



1 . QO 0E + 

!.' . l; ; ■() U (i Oui. -0 1 
•' .'in I) u OIL ■ i. 

2 . 5 E 1 
L . ':■ . OIL- i 
i . 420 57 IE -01 



5 OOOOOOi -0i 

L 01. ■■ i 

L !.! i. Oi ■•■() 1. 

;.,.-: u (i 1. ■■■■ i 

2 U i 'I 



3. 333 333 L ■■■■ i 

3 . LL.L.LL! 0.1 

", L' : *3'* -> 01 
i . 666667E--0 1 
1 . 66 6 66 7 E -0 i 



Results: 



MATRIX IN SYMMETRIC STORAGE MODI.- 



1 . ( i o o o i; : ; •*■ (i o •;;; . q q o (j q (i £ •- q j 
} l 2 5 (i 1. ■■■■ (; 
.L0 1:- ■■■■ C 



"r ■;: -;; "f y~ .... 



LOO 00 it r -0 i 



i 2 U !. 
i 1 . 66666 7 E 



i 
01 

j 1 



• Symmetric to full storage mode 
User entries: 



ORDER OF MATRIX IN FULL STORAGE MODE 



S Y M M E T R I ( :: S T RAG E MATRIX 



1 . 000000 E + 00 

3 . A3-A z - Lii. -Oi. 

2 . 50000 0E-01 
1 .42837 IE- 01 



.! ::. 



'■■: ( 

2 .500 o o o i:: 
2 . o o e; 



•0 i 3 333333!' --Oi 
1 2 i) 01 • 1 
0:1 1 . 666667E -01 



Results: 



MATRIX IN FULL STORAGE MODE 



O0O000E+OO 
50 0E-0 



300 00E-0 
50000 OIL- 
30 00 0E- 
42857 IE 



5. 00 00F- 01 
5. 000 00 0E- 01 
2. 000 0E- 01 
2. 00 00 OF -01 
2. 00 00 OF -01 



3.333333E 
3.333333E 
1 . 66 6 6 67 E 
1.666667E 



■O 1 

• 1 
1 

• 1 
01 



01 
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Chapter O 
Eigen Analysis 



Introduction 

Description 

This section contains several routines for computing the eigenvalues and eigenvectors of real 
variables. Suppose A is a square matrix, and consider the equation Ax = \x, where x is a 
vector and A. is a constant. A number X for which this equation has a non-zero solution vector 
x is called an eigenvalue of the matrix A. The solution x is called an eigenvector of the matrix 
A. 

Routines 

• Eigen - finds all the eigenvalues and eigenvectors of a real general matrix. 

• Compve - finds the complex eigenvector of the real upper-Hessenberg matrix of order n; 
contained in file "Eigen". 

• Realve - finds the real eigenvector of the real upper-Hessenberg matrix in the array A; 
contained in file "Eigen". 

• Scale - scales a general matrix so that the quotient of the absolute sum of the off- 
diagonal elements of column i and the absolute sum of the off-diagonal elements of row i 
lies within certain values; contained in file "Eigen". 

• Symqr - finds the eigenvalues and eigenvectors of a real symmetric matrix. 
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Eigenvalues and Eigenvectors 
of a Real General Matrix 

Description 

This subprogram finds all the eigenvalues and eigenvectors of a real general matrix. The eigen- 
values are computed by the QR double-step method and the eigenvectors by inverse itera- 
tion. 

Program Usage 

Driver Utilization 

• File Name - "EIGEN", disc 1 

The driver "EIGEN" sets up the necessary input parameters for the subprogram Eigen and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Eigen", disc 1 

• Calling Syntax 

CALL Eigen (N, A(*), Evr (*), Evi(*), Vecr(*), Veci(*), Indic(*)) 

• Input Parameters 



N 
A(*) 

Output Parameters 



Order of matrix A. 

Array containing matrix for which eigenvalues and eigen- 
vectors are to be found; dimensioned A(1:N, 1:N). 



A(*) 
Evr(*) 

Evi(*) 

Vecr(*) 

Veci(*) 

Indic(*) 



The original matrix A is destroyed. 

Vector containing the real parts of the n computed eigenval- 
ues; dimensioned Evr(l:N). 

Vector containing the imaginary parts of the n computed 
eigenvalues; dimensioned Evi(l:N). 

Array containing the real components of the normalized 
eigenvector i (for i = 1 to n), corresponding to the eigenval- 
ue stored in Evr(I) and Evi(I); the ith eigenvector is stored in 
column i; dimensioned Vecr(l:N, 1:N). 

Array containing the imaginary components of the normal- 
ized eigenvector i (for i = 1 to n), corresponding to the 
eigenvalue stored in Evr(I) and Evi(I); the ith eigenvector is 
stored in column i; dimensioned Veci(l:N, 1:N). 

Array indicating the success of the subprogram as follows: 

Value of Indic(I) Eigenvalue i Eigenvector i 

not found not found 

1 found not found 

2 found found; 
dimensioned Indic(l:N). 
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• Subprograms Required 

Scale This subprogram scales matrix A so that the quotient of the 

absolute sum of the off-diagonal elements of column i and 
the absojute sum of the off-diagonal elements of row i lies 
within certain bounds. 

Hesqr This subprogram finds all the eigenvalues of a real general 

matrix. 

Realve This subprogram finds the real eigenvector of the real up- 

per-Hessenberg matrix in the array A, corresponding to the 
real eigenvalue stored in Evr(Ivec). The inverse iteration 
method is used. 

Compve This subprogram finds the complex eigenvector of the real 

upper-Hessenberg matrix of order n corresponding to the 
complex eigenvalue with the real part in Evr(Ivec) and the 
corresponding imaginary part in Evi(Ivec). The inverse itera- 
tion method is used in a modified manner to avoid the use 
of complex arithmetic. 

Further explanation of these subprograms is given on pages 5-8 through 5-14. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Eigen." 

N = 

The data may be corrected from the keyboard (e.g., N = 7, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• The Fortran program from which this subprogram has been adapted (Ref. 1), has been 
extensively tested (Ref. 2). Here is a quote from their conclusions: 

"Conclusions. The algorithm is capable of successfully computing 
eigenvalues and eigenvectors of real general matrices even under 
conditions considered unstable. It has the advantage of being com- 
putationally fast, and has the capability of yielding results with as 
much precision as the hardware will permit. The algorithm does not 
break down when presented with a matrix which is not diagonaliza- 
ble; that is, a set of eigenvectors satisfying the eigenequation is com- 
puted regardless of the existence of linearly independent eigenvec- 
tors. However, when a matrix is diagonalizable and degenerate, the 
algorithm does not yield well separated eigenvectors corresponding 
to non-distinct eigenvalues. Another apparent disadvantage is the 
possible indication of completely successful computation (INDIC), 
even in clearly ill-conditioned situations where computational difficul- 
ties are inevitable. This latter property, however, is a common fault 
of other algorithms as well." 
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Methods and Formulae 

This subprogram finds all the eigenvalues and the eigenvectors of a real general matrix of 
order n. 

First, in the subprogram Scale the matrix is scaled so that the corresponding rows and col- 
umns are approximately balanced and then the matrix is normalized so that the value of the 
Euclidian norm of the matrix is equal to one. 

The eigenvalues are computed by the QR double-step method in the subprogram Hesqr. The 
eigenvectors are computed by inverse iteration in the subprogram Realve, for the real eigen- 
values, or in the subprogram Compve, for the complex eigenvalues. 

The elements of the matrix are stored in the two-dimensional array A. The original matrix is 
destroyed by the subprogram. 

Upon output from the subprogram, the real parts of the n computed eigenvalues will be found 
in the first n places of the array Evr and the imaginary parts in the first n places of the array 
Evi. The real components of the normalized eigenvector i (for i = 1, 2, ..., n) corresponding 
to the eigenvalue stored in Evr(I) and Evi(I) will be found in the first n places of the column i 
of the two-dimensional array Veer and the imaginary components in the first n places of the 
column i of the two-dimensional array Veci. 

The real eigenvector is normalized so that the sum of the squares of the components is equal 
to one. The complex eigenvector is normalized so that the component with the largest value 
in modulus has its real part equal to one and the imaginary part equal to zero. 

The array Indie indicates the success of the subprogram Eigen as follows: 

Value of Indic(I) Eigenvalue i Eigenvector i 

not found not found 

1 found not found 

2 found found 

References 

1. Grad, J., and Brebner, M.A., "Algorithm 343: Eigenvalues and Eigenvectors of a Real 
General Matrix." Comm. ACM , 11 (Dec. 1968), pp. 820-826. 

2. Knoble, H.D., "Certification of Algorithm 343: Eigenvalues and Eigenvectors of a Real 
General Matrix." Comm. ACM , 13 (Feb. 1970), pp. 122-124. 

3. Wilkinson, J.H., The Algebraic Eigenvalue Problem, Oxford: Clarendon Press, 1965, 
pp. 86-93. 
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User Instructions 

1. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "EIGEN" and press RUN. 

2. You will be asked to supply entries for the following items: 

- order of real general matrix A 

- elements of matrix A. 

Press CONTINUE after each entry. 

3. The real and imaginary parts of the eigenvalues are then printed as well as the corre- 
sponding real and imaginary parts of the eigenvectors in the following manner: 

a. Vector Evr contains the real parts of the n computed eigenvalues. 

b. Vector Evi contains the corresponding imaginary parts of the n computed eigenval- 
ues. 

c. Matrix Veer contains the real components of the normalized eigenvector i (for i = 1 
to n), corresponding to the eigenvalue stored in Evr(I) and Evi(I); the ith eigenvec- 
tor is stored in column i. 

d. Matrix Veci contains the imaginary components of the normalized eigenvector i (for 
i = 1 to n), corresponding to the eigenvalue stored in Evr(I) and Evi(I); the ith 
eigenvector is stored in column i. 
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Examples 

User entries: 



Results: 



ORDER OF MATRIX A<*>= 3 



MATRIX A 



8.00000 0E + •••• i . E + -S .000 E+ 
4. 0000O0E + 0O 4. 00O00OE + 0O -2. O00000E + 0O 
i .80 0E+0 i -5 . 00E + -7 .00000 0E + 



REAL COMPONENTS OF EIGENVALUES: 
2 . 000OOOE + 0O 2 . 0OO00OE+0O i . 0OO000E+O0 

IMAGINARY COMPONENTS OF THE EIGENVALUES: 
4.0000 E + •■•• 4.00 E +00 0.0 E + 



REAL COMPONENTS OF EIGENVECTORS 
I! CONTAINED IN COLUMNS:!: 

S. OOOOOOE-Oi 5. OOOOOOE-Oi -4 . 082483E -01 
3 . 0S3113E-- 5.6 3 . 053113E- 16 -8 . i 64966E- 01 
1.00 E+ 1.000000 E + -4 . 82483E- 1 



IMAGINARY COMPONENTS OF EIGENVECTORS 
i: CONTAINED IN COLUMNS]. 

5. OOOOOOE-Oi -5. OOOOOOE-Oi O.OOOOOOE+00 
i.OOOOOOE+00 -1 . 000000E+00 O.OOOOOOE+00 
0. OOO0O0E+O0 O.00OO0OE+0O O.OOOOOOE+00 
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User entries: 



Results: 



ORDER OF MATRIX A<*> = 4 



MATRIX A: 

4.000000E+00 ~S. OOOOOOE+00 0. OOOOOOE+00 

3.Q00Q00E+O0 O.OOOOOOE+00 4.000000E+00 

-3.000000E+00 -5.000000E+00 5.000000E+00 

-3.000000E+00 4.000000E+00 O.OOOOOOE+00 

3.000000E+00 O.OOOOOOE+00 5. OOOOOOE+00 

4. 00OO00E+00 



REAL COMPONENTS OF EIGENVALUES: 

i.200000E+Oi 1.000000E+00 i.OOOOOOE+00 
2. OOOOOOE+00 



IMAGINARY COMPONENTS OF THE EIGENVALUES: 

O.OOOOOOE+00 5. OOOOOOE+00 -5. OOOO00E+O 
. 00000 OE+0 



REAL COMPONENTS OF EIGENVECTORS 
[CONTAINED IN COLUMNS.!: 

-5. OOOOOOE-Oi -i.O0O000E+O0 -i.000O00E+00 

5. OOOOOOE-Oi 5. OOOOOOE-Oi 5 . 777 183E-1 6 

5.777i83E-i6 5. 00 0000E-01 -5 . OOOOOOE-Oi 

-2 . 475935E-16 -2 . 475935E-16 -5. OOOOOOE-Oi 

-S. OOOOOOE-Oi i.OOO0O0E+00 i.00OO00E+OO 

5. OOOOOOE-Oi 



:maginary components of EIGENVECTORS 

[CONTAINED IN COLUMNS]: 

. OOOOOOE+00 8 . 253ii8E-i7 -8.253ii8E-i7 

O.000O00E+00 0.0OO0O0E+0O i.OOOOOOE+00 

•i.OOOOOOE+00 O.OOOOOOE+00 0. 0Q0000E+00 

i.OOOOOOE+00 -i.OOOOOOE+00 0. OO0O00E+O0 

0. 00O0O0E+OO O.OOOOOOE+00 O.OOOOOOE+00 
. OOOOOOE+00 
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Subprogram Compve 

This subprogram finds the complex eigenvector of the real upper-Hessenberg matrix of order 
n corresponding to the complex eigenvalue with real part in Evr(Ivec) and imaginary part in 
Evi(Ivec). The inverse iteration method is used in a modified manner to avoid the use of 
complex arithmetic. 

Compve is used in subprogram Eigen which finds the eigenvalues and eigenvectors of a real 
general matrix. 

Subprogram Utilization 

• File Name - contained in file "Eigen", disc 1 

• Calling Syntax 

CALL Compve (N, M, Ivec, A(*), Vecr(*), H(*), Evr(*), Evi(*), Indict*), Subdia(*), 
Work(*), Eps, Ex) 

• Input Parameters 



N 
M 



Ivec 

A(*) 
Vecr(*) 



H(*)> 

Evr(*) 

Evi(*) 

Subdia(*) 

Work(*) 
Eps 

Ex 



Order of matrix A. 

Order of the submatrix obtained by a suitable decomposi- 
tion of the upper-Hessenberg matrix if some subdiagonal 
elements are equal to zero; the value of M is chosen so that 
the last N - M components of the eigenvector are zero. 

Gives the position of the eigenvalues in the arrays Evr and 
Evi for which the corresponding eigenvector is computed. 

Array used for work space. 

Array containing the real components of the normalized 
eigenvector i (for i — 1 to n), corresponding to the eigenval- 
ue stored in Evr(I) and Evi(I); the ith eigenvector is stored in 
column i. 

Array containing upper-Hessenberg matrix from subpro- 
gram Eigen. 

Vector containing the real parts of the eigenvalues. 

Vector containing the imaginary parts of the eigenvalues. 

Vector containing part of upper-Hessenberg matrix from 
subprogram Eigen. 

Array used for work space in inverse iteration process. 

Small positive number that numerically represents zero; 
from subprogram Hesqr. 

2 -39 
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• Output Parameters 

A(*) 
Indic(*) 



H(*)> 



The contents of array A are destroyed. 

Vector indicating the success of the subprogram as follows: 

Value of Indic(I) Eigenvector i 



1 
2 



not found 
found 



Array containing computed eigenvectors; the real parts of 
the first m components of the computed complex eigenvec- 
tor will be found in the first m places of the column whose 
element is Vecr(l, Ivec) and the corresponding imaginary 
parts of the first m components of the complex eigenvector 
will be found in the first m places of the column whose top 
element is Vecr(l, Ivec- 1). 



Methods and Formulae 

This subprogram finds the complex eigenvector of the real upper-Hessenberg matrix of order 
n corresponding to the complex eigenvalue with real part in Evr(Ivec) and imaginary part in 
Evi(Ivec). The inverse iteration method is used in a modified manner to avoid the use of 
complex arithmetic. 

First, a small perturbation of equal eigenvalues is made if necessary, to obtain a full set of 
eigenvectors. Then Gaussian elimination of the upper-Hessenberg matrix 
((H — Fksi*I)*(H — Fksi*I) + (Eta*Eta)*I) in the array A. The row interchanges that occur are 
indicated in the array Iwork. All the multipliers are stored in the first and second subdiagonal 
of array A. 

The inverse iteration is performed on the matrix until the infinite norm of the right-hand side 
vector is greater than the bound defined as 0.01/(N*Ex). Then the residuals are computed 
and the residuals of the two successive steps of the inverse iteration are compared. If the 
infinite norm of the residual vector is greater than the infinite norm of the previous residual 
vector, then the computed eigenvector of the previous step is taken as the final eigenvector. 



1 H(*) is both an input and output parameter. 
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Subprogram Realve 

This subprogram finds the real eigenvector of the real upper-Hessenberg matrix in the array 
A, corresponding to the real eigenvalue stored in Evr(Ivec). The inverse iteration method is 
used. 



Realve is used in subprogram Eigen which finds the eigenvalues and eigenvectors of a real 
general matrix. 

Subprogram Utilization 

• File Name - contained in file "Eigen", disc 1 

• Calling Syntax 

CALL Realve(N, M, Ivec, A(*), Vecr(*), Evr(*), Evi(*), Work(*), Indic(*), Eps, Ex) 

• Input Parameters 



N 
M 



Ivec 

A(*) 

Evr(*) 

Evi(*) 

Work(*) 

Eps 

Ex 



Order of matrix A. 

Order of the submatrix obtained by a suitable decomposi- 
tion of the upper-Hessenberg matrix if some subdiagonal 
elements are equal to zero; the value of M is chosen so that 
the last N - M components of the eigenvector are zero. 

Gives the position of the eigenvalue in the array Evr for 
which the corresponding eigenvector is computed. 

Array containing values for which real eigenvalues of the 
real upper-Hessenberg matrix are to be found. 

Vector containing the real parts of the eigenvalues. 

Vector containing the imaginary parts of the eigenvalues. 

Array used for work space. 

Small positive number that numerically represents zero; 
from subprogram Hesqr. 
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• Output Parameters 
Vecr(*) 

Indic(*) 



Array containing the real components of the normalized 
eigenvector i (for i = 1 to n), corresponding to the eigenval- 
ue stored in Evr(I) and Evi(I); the ith eigenvector is stored in 
column i. 

Vector indicating the success of the subprogram as follows: 
Value of Indic(I) Eigenvector i 



1 
2 



not found 
found 
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Methods and Formulae 

This subprogram finds the real eigenvector of the real upper-Hessenberg matrix in the array 
A, corresponding to the real eigenvalue stored in Evr(Ivec). The inverse iteration method is 
used. 

First, a small perturbation of equal eigenvalues is made if necessary, to obtain a full set of 
eigenvectors. Then Gaussian elimination of the upper-Hessenberg matrix A is employed. All 
row interchanges are indicated in the array Iwork. All the multipliers are stored as the subdi- 
agonal elements of A. 

The inverse iteration is performed on the matrix until the infinite norm of the right-hand side 
vector is greater than the bound defined as 0.01/(N*Ex). Then the residuals are computed 
and the residuals of the two successive steps of the inverse iteration are compared. If the 
infinite norm of the residual vector is greater than the infinite norm of the previous residual 
vector, then the computed eigenvector of the previous step is taken as the final eigenvector. 
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Subprogram Hesqr 

This subprogram finds the eigenvalues of a real general matrix. The original matrix A of order 
n is reduced to upper-Hessenberg form H by means of similarity transformations (Househol- 
der's Method). 

Hesqr is used in subprogram Eigen which finds the eigenvalues and eigenvectors of a real 
general matrix. 

Subprogram Utilization 

• File Name - contained in file "Eigen", disc 1 

• Calling Syntax 

CALL Hesqr(N, A(*), H(*), Evr(*), Evi(*), Subdia(*), Indic(*), Eps, Ex) 

• Input Parameters 



N 
A(*)' 



H(*) 

Ex 



Order of matrix A. 

Array containing general matrix A; if used as part of sub- 
program Eigen, A(*) contains the scaled and normalized 
matrix A outputted from subprogram Scale; dimensioned 
A(1:N, 1:N). 

Array containing original matrix A; dimensioned H(1:N, 
1:N). 

2 -3 9 



• Output Parameters 

A(*)' 
H(*)' 



Evr(*) 

Evi(*) 

Subdia(*) 

Indic(*) 



Eps 



The input array is destroyed. 

Array containing original upper half of matrix H; the special 
vectors used in the definition of the Householder trans- 
formation matrices are stored in the lower part of array H. 

Vector containing the real parts of the n eigenvalues to be 
found; dimensioned Evr(l:N). 

Vector containing the imaginary parts of the n eigenvalues 
to be found; dimensioned Evi(l:N). 

Vector containing parts of input matrix H; dimensioned 

Subdia(l:N). 

Vector indicating the success of the subprogram as follows: 

Value of Indic(I) Eigenvalue i 




1 



not found 
found; 



dimensioned Indic(l:N). 

Small positive number that numerically represents zero in 
the subprogram; Eps = <Euclidian norm of H>*Ex. 



1 A(») and H (*) are both input and output parameters. 
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Methods and Formulae 

This subprogram finds all the eigenvalues of a real general matrix. The original matrix A of 
order n is reduced to the upper-Hessenberg form H by means of similarity transformations 
(Householder Method). The matrix H is preserved in the upper half of the array H and in the 
array Subdia. The special vectors used in the definition of the Householder transformation 
matrices are stored in the lower part of the array H. 

The real parts of the n eigenvalues will be found in the first n places of the array Evr, and the 
imaginary parts in the first n places of the array Evi. The array Indie indicates the success of 
the routine as follows: 

Value of Indic(I) Eigenvalue i 

not found 

1 found 

Eps is a small positive number that numerically represents zero in the program. 
Eps = < Euclidian norm of H>*Ex. 



1 A(t) and H(*) are both input and output parameters. 
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Subprogram Scale 

This subprogram scales a general matrix so that the quotient of the absolute sum of the 
off-diagonal elements of column i and the absolute sum of the off-diagonal elements of row i 
lies within certain values. 

Scale is used in subprogram Eigen which finds the eigenvalues and eigenvectors of a real 
general matrix. 

Subprogram Utilization 

• File Name - contained in file "Eigen", disc 1 

• Calling Syntax 

CALL Scale (IN, A(*), H(*), Prfact(*), Enorm) 

• Input Parameters 

N Order of matrix A. 

A(*)' Array containing matrix to be scaled and normalized so that 

the Euclidian norm is equal to one; dimensioned A(1:N). 

• Output Parameters 

A(*) : Array containing the scaled matrix. 

H(*) Array containing temporary storage for original A(*). 

Prfact(*) Vector containing the scaling factor; the component i of the 

eigenvector obtained by using the scaled matrix must be 
divided by the value found in ith position of Prfact(*). In 
this way, the eigenvector of the non-scaled matrix is 
obtained. 

Enorm Scaling factor; the eigenvalues of the normalized matrix 

must be multiplied by Enorm in order that they become the 
eigenvalues of the non-normalized matrix. 

Methods and Formulae 

This subprogram stores the matrix of the order n from the array A into the array H. Afterward 
the matrix in the array A is scaled so that the quotient of the absolute sum of the off-diagonal 
elements of column i and the absolute sum of the off-diagonal elements of row i lies within the 
values of Boundl and Bound2. The component i of the eigenvector obtained by using the 
scaled matrix must be divided by the value found in Prfact(I) of the array Prfact. In this way 
the eigenvector of the non-scaled matrix is obtained. 

After the matrix is scaled it is normalized so that the value of the Euclidian norm is equal to 
one. If the process of scaling was not successful the original matrix from the array H would be 
stored back into A and the eigenproblem would be solved by using this matrix. The eigenval- 
ues of the normalized matrix must be multiplied by the scalar Enorm in order that they be- 
come the eigenvalues of the non-normalized matrix. 

For more general information, see subprogram Eigen. 

1 A(*) is both an input and output parameter. 
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Eigenvalues and Eigenvectors 
of a Real Symmetric Matrix 

Description 

This subprogram finds the eigenvalues and, upon your request, the eigenvectors of a real 
symmetric matrix. If the matrix is not initially tridiagonal, it is reduced to tridiagonal form by 
Householder's Method. The eigenvalues of the tridiagonal matrix are then calculated by a 
variant of the QR algorithm with origin shifts. 

Program Usage 

Driver Utilization 

• File Name - "SYMQR", disc 1 

The driver "SYMQR" sets up the necessary input parameters for the subprogram Symqr and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Symqr", disc 1 

• Calling Syntax 

CALL Symqr (A(*), D(*), E(*), KO, N, Eps, Abscnv, Vec, Trd) 

• Input Parameters 



A(*) 1 

E(*y 

KO 

N 

Eps 
Abscnv 

Vec 
Trd 



Array containing the following: if the matrix is not initially 
tridiagonal, it is contained in the lower triangle of A(*); if the 
matrix is initially tridiagonal, input from A(*) is not used; 
dimensioned A(1:N, 1:N). 

Vector containing the diagonal elements if the matrix is 
initially tridiagonal; subscripted from 1 to N. 

Vector containing the off-diagonal elements if the matrix is 
initially tridiagonal; subscripted from 1 to N- 1. 

An initial origin shift to be used until the computed shifts 
settle down. 

Order of the matrix, i.e., number of rows or number of col- 
umns in the matrix. 

Convergence tolerance. 

Abscnv = 1 if absolute convergence criterion is to be used. 
Abscnv = if relative convergence criterion is to be used. 
See Special Considerations and Programming Hints for 
further details. 

Vec = 1 if eigenvectors are to be computed. 
Vec = if eigenvectors are not to be computed. 

Trd = 1 if matrix is tridiagonal. 
Trd = if matrix is not tridiagonal. 



1 A(»), D(*), and E(*) are both input and output parameters. 
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• Output Parameters 

M*) 1 If eigenvectors are not requested, the lower triangle of A(*) 

is destroyed while the elements above the diagonal are left 
undisturbed; if eigenvectors are requested, they are re- 
turned in the columns of A(*). 

D(*)' Vector containing the eigenvalues of the matrix. 

E(*) J E(I) contains the number of iterations required to compute 

the approximate eigenvalue D(I). 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Symqr." 
N = Eps = 

You may correct the data from the keyboard (e.g., N = 10, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• The maximum number of iterations allowed per eigenvalue is set at 50. If this is ex- 
ceeded, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Symqr." 

"MAX # OF ITERATIONS EXCEEDED ON EIGENVALUE <n>." 

This is a fatal error. (The program may not be continued from this point, but must be 
restarted.) The program pause allows you to query other variables in the subprogram 
environment. 

• To avoid an excessive number of QR steps, an important consideration when eigenvec- 
tors are computed, the following guidelines should be followed. The convergence toler- 
ance should not be smaller than the data warrants (reference 1, p.xxx). The relative 
convergence criterion should be used only when there are eigenvalues, small compared 
to the elements of the matrix, that are nonetheless determined to high relative accuracy. 

• For best results when there is a wide disparity in the sizes of the elements of the matrix, 
the matrix should be arranged so that the smaller elements appear in the lower right- 
hand corner. 



1 A(*). D( + ), and E(*) are both input and output parameter: 
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Methods and Formulae 

Symqr finds the eigenvalues and, if desired, the eigenvectors of a real symmetric matrix. If the 
matrix is not initially tridiagonal, it is reduced to tridiagonal form by Householder's Method. 
The eigenvalues of the tridiagonal matrix are calculated by a variant of the QR algorithm with 
origin shifts (see reference 2). 

Eigenvectors are calculated by accumulating the products of the transformations used in the 
Householder transformations and the QR steps, a procedure which guarantees a nearly ortho- 
normal set of approximate eigenvectors. 

At each QR step, the eigenvalues of the 2x2 submatrix in the lower right-hand corner are 
computed, and the one nearest the last diagonal element is distinguished. When these num- 
bers settle down, they are used as origin shifts. 

You may choose between absolute and relative convergence criteria. The former accepts the 
last diagonal element as an approximate eigenvalue when the last off-diagonal element be 
small compared to the last two diagonal elements. 

References 

1. Stewart, G.W., "Eigenvalues and Eigenvectors of a Real Symmetric Matrix", Comm. 
ACM (June 1970), pp. 384-386. 

2. Stewart, G.W., "Incorporating Origin Shifts into the Symmetric QR Algorithm for Sym- 
metric Tridiagonal Matrices", Comm. ACM (June 1970), pp. 365-367. 

3. Wilkinson, J.H., The Algebraic Eigenvalue Problem, Oxford: Clarendon Press, 1965. 

User Instructions 

1. Make sure Numerical Analysis flexible disc 1 is inserted correctly into the disc drive. Load 
file "SYMQR" and press RUN. 

2. You will be asked to supply entries for the following items: 

- order of real symmetric matrix A (number of rows or number of columns) 

- elements of matrix A if matrix is not initially tridiagonal 

- elements of vector D if the matrix A is initially tridiagonal where the ith element of D 
is the element (i, i) of matrix A. 

- elements of vector E if the matrix A is initially tridiagonal where the ith element of E 
is the element (i + 1, i) of matrix A. 

- initial origin shift 

- absolute or relative convergence criterion usage 

- convergence tolerance. 

Press CONTINUE after each entry. 

3. The eigenvalues, as well as the eigenvectors, if desired, will be printed. Eigenvector i 
corresponding to eigenvalue i is contained in column i of the eigenvector matrix. 
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Example 

User entries: 



URL I- k Ul- A i *.- ;: 4 



MATRIX A : 

':•('! U I.) f. + IJ 

i. . oil ( i 1 ■*• U (j 
:i. . 0OOOOE < 00 
i fi i.i U (I Oi Ml (i 
:'i. . (5 E + 
4 0E-* 00 



4.00 E '<■ 
i Q (i i) i * 
:t. . E ■*■ (i 
4 E ■*• 
i (iQ (i l.i Hi +0U 



1 U !" + 

'::. o (i n o :) or '• o o 

1 . i) I:; (• U 
:! f i (i .) \ :■■ 
2.00 E + 



IN II IAi.. ORIGIN SHIFT-- 
LUNUERGENGE TOLERANCE- i. . IE S 



Results: 



EIGENVALUES 



:t 00000 E+Oi 5 . 000 (JOE + 00 

2 OOO 0U0E-+00 



i . OOOOOOEHiO 



e I G t N v 1 E C T R S : i: C N. T a I n E I) I N C 1... M N s :s 



324SSSE"0:j 
0000 00E +00 
0710 68E -0 i 
32 4 S3 3 E-- 01 
:i.6"-'273i-. - Oi. 
0710 68E -Oi 



3. i 62 278 E- 01 

6 324555E--01 

3 . 833642 E- 17 

■2 53922SE---16 

-6. 32433SE--01 



-7.071068E-- 01 

3 . 162 2 7 BE -01 

3 162278E-01 

-7. 071068E-01 

■3.S3 0738E-16 
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Chapter 6 
Interpolation 



Introduction 

Description 

This section has three interpolation programs. 

Routines 

• Dvdfc - computes a table of confluent divided differences. 

• Bnewt - performs Newton interpolation with backward divided differences; used in con- 
junction with Dvdfc. 

• Fnewt - performs a Newton interpolation with forward divided differences; used in con- 
junction with Dvdfc. 

• Spline - computes a curve s(x) that passes through n data points (x i5 y t ). 

• Cheby - fits the tabular function Y(X) (given as m points (X, Y)) by a polynomial 

n 

P = 2 AjX' 

i = 
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Confluent Divided Differences 

Description 

Given a set of real numbers, {Xj, X 2 , ..., X k }, and a corresponding set of function values, {V^, 
^2, •••, V k }, the forward divided differences are defined as follows: 

Zero order differences are 

Af(XJ = V n , n = 1, ..., k 

First order differences are 

Af(X n , X n + 1 ) = V n + 1 - V n )/(X n + 1 - XJ. n = 1, ..., k - 1. 

Higher order differences are defined in terms of lower order differences: 
Af(X b X i + 1 , ..., X n ) = (Af(X i + 1 , ..., XJ- Af(X h ..., X^))/^ - X,), n = i + 2, ..., k. 

The differences may be displayed in the form of a table: 



Xj 


V: 










Af(X l5 


X 2 ) 


x 2 


V 2 




Af(X : , X 2 , X 3 ) 






Af(X 2 , 


X 3 ) Af(X lt X 2 , X 3 , X 4 


x 3 


v 3 




Af(X 2 , X 3 ,X 4 ) 






Af(X 3 , 


X 4 ) 


x 4 


v 4 






For example: 








X 


f(X) 






2.0 


1.11 


0.15 




2.2 


1.14 




3.33E - 2 






0.17 


- 8.08E - 2 


2.5 


1.19 


0.12 


- 5.56E - 2 


3.1 


1.26 







This program calculates f(Xx, X 2 , ..., X n ) for any integral value of n in the interval [2, k]. 

Program Usage 

Driver Utilization 

• File Name - "DVDFC", disc 2 

The driver "DVDFC" sets up the necessary input parameters for the subprogram Dvdfc and 
prints out the resulting outputs. Either subprogram Bnewt, a Newton interpolator with back- 
ward divided differences, or subprogram Fnewt, a Newton interpolator with forward divided 
differences, may be called to locate points of interpolation. 
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Subprogiram Utilization 

• File Name - "Dvdfc", disc 2 

• Calling Syntax 

FN Dvdfc (N, X(*), V(*), B(*)) 

• Input Parameters 



N 
X(*) 



V(*) 
B(*) 1 



• Output Parameters 

B(*y 

FN Dvdfc 



The number of data points to be used in the calculation. 

Vector containing the initial x-values; the values x s need not 
be distinct or in any special order, but once the vector X is 
chosen it will determine the interpretation of B(*) and V(*); 
subscripted from 1 to at least N. 

Vector containing the values of the functions f(X); sub- 
scripted from 1 to at least N. 

Vector containing backward differences. When N — 1, the 
state of B(*) is irrelevant. When N is greater than 1, B; must 
contain Af(X b ..., X n _i) for i = 1, 2, ..., n - 1 before Dvdfc 
is called. 



Vector containing backward divided differences. After Dvdfc 
is called, you will find Bj - Af(X,, X i + 1 , ..., X n _!, X n ) for 
i = 1,2, ..., n - 1, n. 

Value of the forward divided difference Af(X l5 X 2 , ..., X n ). 



Special Considerations and Programming Hints 

• The values X need not be distinct or in any special order, but once the array X is chosen, 
it will fix the interpretation of B(*) and V(*). If X 1; X 2 , ..., X n are in monotonic order, 
then the effect of roundoff upon any nth divided difference is not more than would be 
caused by perturbing each f(X s )) by n units at most in its last significant place. But if the 
X's are not in monotonic order, the error can be catastrophic if some of the divided 
differences are relatively large. 

• The following program segment is an example of how Dvdfc can be used to construct a 
table of forward or backward differences: 

FOR I = 1 TO N 

X(I) = 
V(I) = 

F(I) = FNDvdfc(I, X(*), V(*), B(*)) 
NEXT I 

The array F can be used in subprogram Fnewt or the array B in Bnewt which are also 
contained in the NUMERICAL ANALYSIS package. 



1 BU) is both an input and output parameter. 
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Methods and Formulae 

A full explanation of the algorithm employed may be found in reference 1. 

References 

1. Kahan, W. and Farras, I., "Algorithm 167: Calculation of Confluent Divided Differ- 
ences", Comm. ACM (April 1963). 

2. Thacher, H., "Certification of Algorithm 167: Calculation of Confluent Divided Differ- 
ences", Collected Algorithms from ACM, pp. 167-168. 

User Instructions 

1. Make sure Numerical Analysis flexible disc 2 is inserted correctly into the disc drive. Load 
file "DVDFC" and press RUN. 

2. You will be asked to supply entries for the following items: 

- number of data points to be used in the calculation 

- initial x-value for each data point (matrix X) 

- initial value of f(x) for each data point (matrix V) 

- choice of backward or forward divided differences 

- domain value (z) of the point to be interpolated. 
Press CONTINUE after each entry. 

3. The value of the interpolated point and its derivative, as well as an error estimate is 
printed. 
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Subprogram Bnewt 

This subprogram performs Newton interpolation with backward divided differences. Bnewt 
may be used in conjunction with subprogram Dvdfc, a routine for calculating confluent di- 
vided differences. 

Subprogram Utilization 

• File Name - "Bnewt", disc 2 

• Calling Syntax 

CALL Bnewt (Z, N, X(*), B(*), P, D, E) 

• Input Parameters 

Z Domain value of the point to be interpolated. 

N Number of data points used. 

X(*) Vector containing the initial x-values; the values Xj need not 

be distinct nor in any special order, but the components 
must correspond to the values in B(*); subscripted from 1 to 
at least N. 

B(*) Vector containing the backward divided differences 

Bj = Af(X i; X i + 1 , ..., X n ) for i = 1 to n; subscripted from 1 
to at least N. 

• Output Parameters 

P Value of the following polynomial in Z of degree N - 1 at 

most: 

B(N) + (Z - X(N))*{B(N - 1) + (Z - X(N - 1){B(N - 2)+ ... +(Z - X(2))B(1)}...} 

This polynomial is an interpolation polynomial which 
would, but for rounding errors, match the values of the 
function f(X) and any of its derivatives that subprogram 
Dvdfc might have been given. 

D Value of the derivative of P. 

E Estimated maximum error in P caused by roundoff during 

the execution of Bnewt. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Bnewt." 

N = 
You may correct the data from the keyboard (e.g., N = 20, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

Methods and Formulae 

A full explanation of the algorithm employed may be found in reference 1, page 6-4. 
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Subprogram Fnewt 

This subprogram performs a Newton interpolation with forward divided differences. Fnewt 
may be used in conjunction with subprogram Dvdfc, a routine for calculating confluent di- 
vided differences. 

Subprogram Utilization 

• File Name - "Fnewt", disc 2 

• Calling Syntax 

CALL Fnewt (Z, N, X(*), F(*), R, D, E) 

• Input Parameters 

Z Domain value of the point to be interpolated. 

N Number of data points used. 

X(*) Vector containing the initial x-values; the values X s need not 

be distinct nor in. any special order, but the components 
must correspond to the values in B(*); subscripted from 1 to 
at least N. 

F(*) Vector containing the forward divided differences 

Fj = Af(Xj, X i + 1 , ..., X n ) for i = 1 to n; subscripted from 1 
to at least N. 

• Output Parameters 

R Value of the following polynomial in Z of degree N - 1 at 

most: 

F(1) + (Z - X(1))*{F(2) + (Z - X(2)){F(3)+ ... + (Z - X(N - 1))F(N)}...} 

This polynomial is an interpolation polynomial which 
would, except for rounding errors, match the values of the 
function f(X) and any of its derivatives that subprogram 
Dvdfc might have been given. 

D Value of the derivative of R. 

E Estimated maximum error in R caused by roundoff during 

the execution of Fnewt. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Fnewt." 
N = 

You may correct the data from the keyboard (e.g., N = 5, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 
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Methods and Formulae 

A full explanation of the algorithm employed may be found in reference 1, page 6-4. 



Examples 

User entries: 



NUMBER OF DATA POINTS= ii 

X ( * ) V ( * :> 



Results: 



i 


••■•5 


000000E+00 


-5 


0OO0O0E+OO 




•-3 


OOOOOOE+0 


-3 


OOOOOOE+00 


-.':< 
.<}. 


•••• 1 
1 


OOOOOOE + 0- 

000 0E +00 


•••• 1 
i 


0OO0O0E+0O 

OOOOOOE+00 


'.'.'.' 


3 


0O00OOE+0O 


3 


OOOOOOE+0 


{..* 


c; 


OOOOOOE+0 




OOOOOOE+0 


,-, 


*"/' 


O0 00OOE+0 


i»i 


OOOOOOE+0 


& 


9 


OOOOOOE+0 


9 


OOOOOOE+00 


o 


i 


iO0OOOE+0i 


i 


10 OE+Oi 


i. 


i 


30 00 0E+01 


l 


30 0GOOE+01 


1 i 


i 


50000 0E+ 01 


l 


500000E+01 



B<# 



*> 



1 I' 



•5 . OOOOOOE + 00 
3 . OOOOOOE + 
■1 . OOOOOOE+00 
i . OOOOOOE + 00 
3 . E + 
'::■ . O0OO00E + 00 
7 . OOOOOOE+00 
9 . OOOOOOE+00 

i :i. o o o o o e + o .1 

1 30 000E+01 
1 S00 0GE+0I 



■S. 0OOO00E+O0 
i . O0O00OE + O0 


(j 





fl 



E + 
OOOOOOE+0 
E + 
000 00 0E+0 
E + 
00O0O0E+0O 
E + 
000 0E+0 
E + 



NEWTON INTERPOLATION WITH 

! ;: ' O R W A R D I) I U I D E I) D I F l ;: ' E R E N O E S ■■ 



I N T E R P O !... A T E D P T . - 2.00 E + 
I) E R I 7 A T I 7 E. ■■■■■ .1. . () E. + 
ERRORS 3 450 000E--07 

I! fl "f ' E R P O L A "f E D ! :! T . ■■■■■■■■ ■■■■ 2 E + 

i) i::: r :r y a t j: y e ■■■■■■ i. o o o o o o e + o o 

ERROR- i . 630 0E- 07 

I! N "!' p: R P O i .. h T e v> p t . ■■■■■■■■ X . 2 E *■ L i 
0\:\i \ ' ",'"' ; j ■.,'!' :;: 1 ii if :. •'■ fi 
ERR OP -■■■ 7 . 430 0E/-0 7 
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You enter the value of Z, Fnewt computes the rest: 



User entrie; 



NUMBER OF DATA POINTS 



13 



X ( $ > 



v 1 (* 



i 


. 00O000lii>O0 





< .. 


i . o o o o o o f * o o 


7 


3 

■ ■." 


7 t.)v +0 
3.00 E v- ('; 
4.00 i'J E '< 
5 . OO00 0E+0O 

f ) (1 i) (1 |i : -f (.1 () 


R 
1. 


8 


7 (J 00 OH 01::.+ (J 


9 


9 


8 . () Hi! + (5 


i 'i 


1 


9 . O00OOOE+OO 


''.' 


i. 1 


1 . O0000OE+01 


1'"' 
"> 


>. 2 


i i (I n F + (i 1 


;■> 


i. 3 


1 200 OF + 01 






Results: 



B<*> 



f i !. ■''■ 
SB 919 0E -01 
() 01 - 1 
071068E- 01 
660254E-01 
65 9 27 BE -01 
1 •*■ 
67927BE-01 
66 254E-01 
071068E-01 
00O0OOE--01 
58819 017-01 
li!> 

I" ( * ) 



I; Oi *■ 



10 
1 .1 
1 2 
13 



■\'.'f 


5881 9 0E 


■-• 1 


'.'.) 


. E 


•- 1 


'/' 


0710 68F 


- i 


p 


6 6 (.' 2 b 4 1::. 


•■•• i 


9 


659 2 5 BE 


-• i 


1 


. n fl 1"' 


+ 


9 


659? 5 BE 


-■ 1 


!'.) 


6602S4F 


■01 


''.'■' 


0710 68 F 


■■■ i 


C 


E 


■••• i 


('... 


588 19 OF 


-Oi 





iii. 


+ 






(.1 (J 1::. 


* 


'■■' 


58B190E 


- 1 


3 


8 190 0E 


- 3 




739367E 


-03 


9 


675833E 


- 5 


B 


1 5 liii 


- 6 


-..'> 


10 8333 Iii: 


- 7 


1 


liii 


- B 


4 


Oi 150 BE 


• i 


1 


102293E 


- 1 2 


1 


488 09SE 


- 1 2 


3 


8S8 025E 


• • 1 3 


6 


43 411;:: 


••■ 1 4 



n iii: w t o n :i: n t iii: r p o i... a t i o n w i t h 

F O R W A R D D I V I D liii D V> 1 F F liii R liii N C Iii. 8 



Z™ 1 . 5 

I N T liii R P O I... A T E I) P T . -= 3 . B 2 6 8 3 5 F - 1 
D liii R I U A T I V E ■■■■■■■ 2.418710 ii: - i 
ERRORS 2. 32611 BE- 08 



/.- 11 . 2 b 

I N T Iii. R P O i... A T Ei D P T . ~= 1.950901 E •■•• : 

i) i::: r i m a t i v f - • - 2 . b 6 ? 6 9 9 f •••• 1 

ERROR- 2. B939B4E-07 
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Cubic Spline Interpolator 

Description 

This subprogram computes a curve s(x) that passes through n data points (x s , y^ that you 
supply. 

The integral /**" s(t)dt is calculated, as well as s(x) and s'(x) for any point x in the interval [x x , 
x n ]- 

Program Usage 

Driver Utilization 

• File Name - "SPLINE", disc 2 

The driver "SPLINE" sets up the necessary input parameters for the subprogram Spline and 
prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Spline", disc 2 

• Calling Syntax 

CALL Spline (N, Narg, X(*), Y(*), Domain (*), Func(*), Deriv(*), Int, Eps) 

• Input Parameters 



N 
Narg 

X(*) 

Y(*) 
Domain (*) 

Eps 

• Output Parameters 
Func(*) 

Deriv(*) 

Int 



Number of data points. 

Number of arguments for which the derivative or functional 
value is to be computed; Narg = means that only the in- 
tegral is to be calculated. 

Vector containing domain values of the data points; sub- 
scripted from 1 to N; values should be entered in increasing 
order. 

Vector containing range values of the data points; sub- 
scripted from 1 to N. 

Vector containing the domain values for which the deriva- 
tive or functional value is to be computed; subscripted from 
1 to Narg; arguments may be entered in any order, but 
these values must be contained in the interval [x lf x n ]. 

Error tolerance in iterative solution of simultaneous equa- 
tions. 



Vector containing the interpolated function values for out- 
put arguments in Domain(*); subscripted from 1 to Narg. 

Vector containing the derivative values for output argu- 
ments in Domain (*); subscripted from 1 to Narg. 

Integral /x n s(x)dx. 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Spline." 

Eps = N = 

The data may be corrected from the keyboard (e.g., Eps = IE - 6, EXECUTE). When 
CONTINUE is pressed, the program will resume execution at the next line. 

• The arguments for which the derivative or functional values are to be computed must lie 
in the interval [x 1; x n ]. If an argument is found outside this range, the following error 
message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Spline." 

"ARGUMENT OUT OF BOUNDS." 

X_(l) = X_(N) = Domain (<i>) = 

The program may be saved in the following way: 

a. Type: Domain (<i>) = < a permissible value > 

b. Press: EXECUTE 

c. Type: CONT Corrector (Type each letter. Do not use the CONTINUE key. ) 

d. Press: EXECUTE 

• The data points (x b yj, for i = 1 to n should be entered in increasing order of the x h 
where the Xj are discrete and x ( < x i + 1 for i = 1 to n - 1. The output arguments tj, for 
j = 1 to Narg may be entered in any order. 

• The error factor of the solutions is approximately equal to h 4 for the integral, h 3 for the 
functional values and h 2 for the derivative, where h is the average interval size. 

• X(*), Y(*), Domain(*), Deriv(*) and Func(*) must be dimensioned in the calling 
program. 
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Methods and Formulae 

This subprogram computes a curve s(x) that passes through the n data points (x b y,) that you 
supply and computes certain information at any point t ; on the curve as long as tj is in the 
interval [x t , x n ]. The information that can be computed is the integral over the interval and 
the derivative or functional value at any point on the interval. 

The method implemented fits a curve through the points and integrates, differentiates and 
interpolates that curve. The curve used is the cubic natural spline which derives its name from 
a draftsman's mechanical spline. If the spline is considered as a function represented by s(x), 
then the second derivative s"(x) approximates the curvature. For the curve through data 
points (x b yj), (x 2 , y 2 ), ..., (x n , y n ) we want jC" (s"(x)) 2 dx to be minimized in order to 
achieve the "smoothest" curve. 

The spline function with minimum curvature has cubic polynomials between adjacent data 
points. Adjacent polynomials are joined continuously with continuous first and second deriva- 
tives as well as s"(xi) = s"(x n ) = 0. 

The procedure to determine s(x) involves the iterative solution of a set of simultaneous linear 
equations by Young's method of successive overrelaxation. You can specify the accuracy to 
which these equations are solved. 

The formulae employed are: 

jT ,+l s(x)dx= | ~2 (x i+i ~ x i) ( yi + yi+i) - 24 (Xi+i ~ Xi)3 ^ s " (Xi) + s "( x i+i^J 

s(tj) = y, + (tj - x,) ( * + 1 Z* ) + (^ - x,) (t ; - x i + 1 ) \ [s"(x,) + s"(x i + 1 ) + s"(tj)] 

s'(tj) = v Vi + 1 I V ' + \ [(tj - X.) + (tj - X i + 1 )] (S"(X,) + S"(tj)) + T 
x i + 1 x i ° u 



/ s"(x 1 + 1 ) - s"(x.) \ 

( t j- Xi )(tj-x i+1 ) ^ — x^t^^ ) 
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References 

1. Ralston and Wilf, Mathematical Methods for Digital Computers, Vol. II, New York: John 
Wiley and Sons, 1967, pp. 156-158. 

2. Greville, T.N.E., Ed., "Proceedings of An Advanced Seminar Conducted by the 
Mathematics Research Center", U.S. Army, University of Wisconsin, Madison, October 
7-9, 1968, Theory and Applications of Spline Functions, New York, London: Academic 
Press, 1969, pp. 156-167. 

User Instructions 

1. Make sure Numerical Analysis flexible disc 2 is inserted correctly into the disc drive. Load 
file "SPLINE" and press RUN. 

2. You will be asked to supply entries for the following items: 

- number of data points to be used in the calculation 

- number of domain points or number of arguments for which the derivative or in- 
terpolated functional value is to be computed. If only the value of the integral is de- 
sired, enter 0. 

- error tolerance for the iterative solution of the simultaneous equations 

- x-value (in increasing order) and y-value for each data point 

- domain values to be evaluated in the cubic spline. 
Press CONTINUE after each entry. 

3. The program will print the resulting interpolated functional values, the values of the 
derivative and the value of the integral over the interval [Xj, x n ]. 
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Example 

User entries: 



# OF DATA POINTS^ S 

# OF DOMAIN POINTS'* 4 
ERROR TOLERANCE* i . E-6 



DATA POINTS 



X( i )»■••• 3 . OOOOOOE+00 Y< 1 ) 

X( 2) ~-i. OOOOOOE+00 Y< 2> 

X< 3)~ 2. OOOOOOE+00 Y< 3) 

X< 4 )•.::: 4. OOOOOOE + 00 Y( 4) 

X( S> :: " 7 OOOOOOE + 00 Y< 5) 






o 




+ t 
+ I 
+ 0! 
+ o 
+ 



DOHA IN VALUES: 

D o m a :i n (. i > :::: - 2 . S E + 

Donain( 2)= 0. OOOOOOE+00 

Do«a:i.n( 3 > == 2 . SOOO 0E + 00 

Dc.Ma:i.n< 4 > "= % .OOOOOOE + 00 



Results: 



: N T E G R A L F R h •■•■ 3.00 E + T O 
7.00 000 E + ~ 4 . 9 8 4962 E + : 



DOMAIN VALUES 



FUNCTION VALUES 



2 ''. fj iin IF + 
. 00 0F.+ 
2.5000 E + 
'•■; . 00OO0OE+0O 



3.69i8iOE+00 
7 . 7S223SE+0 
6.67 2S9E+0 
1 . 3435S0F+0 



DERIVATIVE VALUES 

i . 399138E+00 

i . 373946E+0 

••-3 . 141092E+0 

S. 306SiE-0i 
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Chebyschev Polynomial 

Description 

This subprogram fits the tabular function Y(X) (given as m points (X, Y)) by a polynomial 

n 

P = S AsX 1 . 

i - 

This polynomial is the best polynomial approximation of Y(X) in the Chebyschev sense. 

Program Usage 

Driver Utilization 

• File Name - "CHEBY", disc 2 

The driver "CHEBY" sets up the necessary input parameters for the subprogram Cheby and 
prints out the coefficients of the resulting polynomial. 

Subprogram Utilization 

• File Name 

• Calling Syntax 

CALL Cheby (M, N, X(*), Y(*), A()) 

• Input Parameters 

M Number of data points to be used in the computation. 

N Degree of the approximating polynomial; for valid results, 

M> (N + 1). 

X(*) Vector containing the x-values of the data points; sub- 

scripted from 1 to M; the Xj must be entered in increasing 
order. 

Y(*) Vector containing the y-values of the data points; sub- 

scripted from 1 to M. 

• Output Parameters 

A(*) Vector containing the coefficients of the polynomial approx- 

imation; subscripted from to N; 
P(X) = A + A X X + A 2 X 2 +...+ A n X n 
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Special Considerations and Programming Hints 

• Although this procedure is an implementation of a finite algorithm, roundoff errors may 
give rise to cyclic changes of the reference set causing the procedure to fail to terminate. 

• X(*), Y(*) and A(*) must be dimensioned in the calling program: DIM X(1:M), Y(1:M), 
A(0:N) 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Cheby." 

"# OF DATA POINTS MUST BE GREATER THAN DEG. OF POLY. + 1." 

M = N = 

The data may be corrected from the keyboard (e.g., M = 6, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

Methods and Formulae 

See the two references for a complete description of the methods and formulae employed in 
subprogram Cheby. 

References 

1. Newhouse, Albert, "Chebyschev Curve Fit", Comm. ACM 5 (May 1962), p. 281. 

2. Stiefel, E., Numerical Methods of Tchebysheff Approximation, U. of Wisconsin Press, 
1959, pp. 217-232. 

User Instructions 

1. Make sure Numerical Analysis flexible disc 2 is inserted correctly into the disc drive. Load 
file "CHEBY" and press RUN. 

2. You will be asked to supply entries for the following items: 

- number of data points to be used in the calculation 

- degree of the polynomial desired for the approximation (for reasonable results, 
(< degree of polynomial) < (< number of data points > - 1) 

— x and y values for each data point. 
Press CONTINUE after each entry. 

3. The coefficients of the Chebyschev Polynomial will be printed where 
P(X) = A + AjX + A 2 X 2 +...+ A n X n . 
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Examples 

User entries: 



Results: 



# Ox DATA POT Mi's = 



DECREE GE POLYNOMIAL APPROXIMATION^ 



DATA POINTS: 

X ( i ) "" • •• S . il E ■*■ Y •. :i ) ■■■■ 2 . 5 E ■* 

X ( 2 > ~" S . (J E •♦• Y ( 2 > ™ 2 ':'■■ n () () E * 

X ( 3 ) "" ••- 4 . i.) E + Y ( 3 ) " :: 1 . 6 E * 

X( 4)- ::: 4 . OOOOOOE-f 00 Y( 4)-- i . c>(.i 00 0E + 

X( S) ■="•■•■ 3 000 OEM) Y( 5 > ;::; 9.0000001':::* 

X< 6)"" 3 . 00 0E + Y( 6 ) :: » 9 OF* 

X ( 7 ) •■=" ■•■■2 . E + Y \ ? ) ■■- 4.00 E * 

X ( 8 ) ■■■■ 2.00 o i;;: t 00 Y ( 8 ) ~= 4 E ■* 

COFFFICTENTB OF GHEBYSCHE'J POLYNOMIAL : 

i. WIT I.' Pi P ( : ' i = :: Ai. ) i [•,•; .1. ) *X !-A ( E)*X '-?f A ( 3 J*X '" Lf 

A < ) - :: . E •* 
A ( i ) ■■■■ 0.00 E *• 
A ( 2> ■■■■ i . E + 
A ( 3 ) "- . E •*• 



01 
(Li 
():'i 








User entries: 



OF DATA POINTS 



Results: 



DEGREE OF PGLYNOhl Ai... APPROXIMATION- 



DATA POINTS: 

X( i )■■■■■■■ ■■ 4 . OOE-MJO Y( i. ) ™ •-■<':• , OEM) 1. 

X ( 2 ) ■■■■■ ■- 2.0000 E + 00 Y i 2 ) -= •■■■ 8 . () E ■*■ 

X( 3>"---i . 00 000 0E*- 00 Y( 3)== :i. . 23560 OEM) 

X( 4)":: 3. . OEM) Y< A)~ i . 00 OEM) 

X( S )■■■■ 3. OO000OEMH3 Y( 5 X ::: 2 . SO 00 OEM) 1 

X( 6)"- S. 00 E+ Y( 6)~ i . 2 0E + 02 



co e f f i c :i: e: n t s o f c h e b y s c h e v p o e y n o m i a l = 

i: WHERE P ( X ) :::: A < > ••:■ A ( .1. ) *X + A ( 2 > *X A 2->- A < 3 ) *X A 3+ . 

AC 0)=-3. 071429E-01 

A( i)---i. 35714300 i 

A( 2)- 3.571429E-02 

A< 3)=- 9.571429E -Oi 
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Chapter 7 
Functions 



Introduction 

Description 

This section contains several mathematical function subprograms that may be used occa- 
sionally. Since the use of these functions may be frequent and well-embedded in a main 
program, no drivers have been supplied with these subprograms. Their usage is explained in 
the following pages. You must supply the calling program. 

Reference 1 below was used as the "cookbook" for these functions. Reference 2 may also be 
useful as a programming reference. 

Routines 

• FN Cosh - hyperbolic cosine 

• FN Sinh - hyperbolic sine 

• FN Tanh - hyperbolic tangent 

• FN Gamma - gamma function 

• FN Lgamma - log gamma function 

• Cadd - addition of complex numbers 

• Cmult - multiplication of complex numbers 

• Cdivid - division of complex numbers 

• Csqrt - square root of a complex number 

• Cexp - exponential value of a complex number 

• Clog - natural logarithm of a complex number 

• Cabs - absolute value of a complex number 

• Cinv - inverse of a complex number 

• Ccos - cosine of a complex number 

• Csin - sine of a complex number 

• Ctan - Tangent of a complex number 

• Ccosh - Hyperbolic cosine of a complex number 

• Csinh - Hyperbolic sine of a complex number 

• Ctanh - Hyperbolic tangent of a complex number 

• Rec_pol - rectangular to polar conversion 

• PoLrec - polar to rectangular conversion 

• Polyev - evaluation of a complex polynomial 
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References 

1. Hart, J., et.al. Computer Approximations, New York: John Wiley and Sons, Inc., 1968. 

2. Cody, William J., Jr. and Waite, William. Software Manual for the Elementary Func- 
tions, Englewood Cliffs, N.J.: Prentice Hall, Inc., 1980. 



Hyperbolic Cosine 

Description 

This subprogram calculates the value of the hyperbolic cosine function. 

Program Usage 

Subprogram Utilization 

• File Name - "Cosh", disc 2 

• Calling Syntax 
FN Cosh(X) 

• Input Parameter 



X 
• Output Parameter 
FN Cosh(X) 

Method and Formula 

The standard formula used is: 



Domain value of the function. 



Value of the hyperbolic cosine at x. 



Cosh(x) = (e x + e~ x )/2 
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Hyperbolic Sine 

Description 

This subprogram calculates the value of the hyperbolic sine function. 

Program Usage 

Subprogram Utilization 

• File Name - "Sinh", disc 2 

• Calling Syntax 
FN Sinh(X) 

• Input Parameter 

X Domain value of the function 

• Output Parameter 

FN Sinh(X) Value of hyperbolic sine at x. 

Methods and Formulae 

The hyperbolic sine function may be defined as follows: 

Sinh(x) = (e x - e~ x )/2 

However, this formula is not appropriate for a computer approximation. Instead, the positive 
domain has been segmented into four sections: 

1. x «s 0.5: Here significance would be lost in the subtraction if the above standard formu- 
la was used. Instead, the following polynomial approximation is employed: 

Sinh(x) - .008420538263*x 5 + .166656747807*x 3 + 1.00000034157 

2. .5 < x ss 1: Sinh(x) - |[D(x) + D(x)/(D(x) + 1)] where D(x) = e x - 1 

3. 1 < x «= 10: Here, the standard formula is used. 

4. x > 10: Here, e" x has lost all significance compared with e x , so we use: Sinh(x) ~ e x /2 

The identity Sinh( -x) = - Sinh(x) is used for negative domain values. 
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Hyperbolic Tangent 

Description 

This subprogram calculates the value of the hyperbolic tangent function. 

Program Usage 

Subprogram Utilization 

• File Name- "Tanh", disc 2 

• Calling Syntax 
FN Tanh(X) 

• Input Parameter 

X Domain value of the function. 

• Output Parameter 

FN Tanh(X) Value of the hyperbolic tangent function. 

Methods and Formulae 

The hyperbolic tangent function may be defined as follows: 

Tanh(x) = (e x - e~ x )/(e x + e~ x ) 

But this formula, by itself, is not appropriate for a computer approximation. 

Instead, the positive domain has been segmented into three sections: 

1. |x| < 1: Here significance would be lost in the subtraction if the above standard formula 
was used, so Tanh(x) « D(2x)/(2 + D(2x)) where D(x) = e x - 1. 

2. 1 =s |x| < 10: Here, the standard formula is used. 

3. |x| > 10: Here, e x has lost all significance compared with e x , so Tanh(x) ~ 1. 

The identity Tanh( - x) = - Tanh(x) allows for the absolute values used above. 
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Gamma Function 

Description 

For a given argument x, (x = 0, - 1, - 2, ...) function subprogram FN Gamma computes 
the value of the gamma function. For an argument x > 70.957, the log gamma function, FN 
Lgamma, should be used to avoid machine overflow. 

Program Usage 

Subprogram Utilization 

• File Name - "Gamma", disc 2 

• Calling Syntax - FN Gamma(X) 

• Input Parameter 

X Function argument; X =£0, - 1, - 2, - 3, ...; arguments 

must be in the range [-69, 70.957] to avoid machine over- 
flow. 

• Output Parameter 

FN Gamma(X) Value of the gamma function evaluated at x. 

Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects an 
argument equal to zero or a negative integer, the following error message is printed and 
the program pauses: 

"ERROR IN FUNCTION Gamma." 

"ARGUMENT VALUE IS EITHER OR A NEGATIVE INTEGER." 

Unless you are willing to change your argument value, there is no recovery from this 
error since the function is not able to handle zero or negative integer arguments. 

• If the subprogram detects an argument which will cause machine overflow, the following 
error message is printed and the program pauses: 

"ERROR IN FUNCTION Gamma." 
"ARGUMENT VALUE OUT OF RANGE." 
ARGUMENT = 

Again, unless you are willing to change your argument value, there is no recovery from 
this error. The argument must be between - 69 and 70.957. 
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Methods and Formulae 

For a given argument x, (x =£0, - 1 ? - 2, ...) this function subprogram computes the value 

n-l 

of the gamma function of x, F(x). The recurrence T(Z + n) = n (Z + k)r(Z) allows the 

k = 

computation of T(x) to be reduced to the computation of T(2 + x), with =s x «= 1. 2 is 
chosen because of the poles at zero and the negative integers. 

Then 

n -1 

T(x) = n (x + 2 + k)r(x + 2) for n > 0, 

k = 



r < x + 2 > for n < 0. 



n (x + 2 - k) 

k=l 



For large x (in absolute value), the above computation is time consuming, and it is more 
economical to use the log gamma function. 
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Log Gamma Function 

Description 

For a given real, non-negative argument x, function subprogram FN Lgamma computes the 
value of the natural logarithm of the absolute value of the gamma function. 

Program Usage 

Subprogram Utilization 

• File Name - "Lgamma", disc 2 

• Calling Syntax 
FN Lgamma(X) 

• Input Parameter 

X Domain value of the function; X must be a real, non- 

negative number. 

• Output Parameter 

FN Lgamma(X) Value of the log gamma function evaluated at x. 

Method and Formula 

For a given real, non-negative x, this function subprogram computes the value of the log 
gamma function of x. 

We use the Stirling form: 

LOG(r(x)) = (x - Vfe)LOG(x) - x + LOG( V2^J + <|>(x) 
where c|>(x) is a rational approximation of the form: 

<MX) = ^Rn, m (l/X 2 ) 
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Addition with Complex Numbers 

Description 

Subprogram Cadd adds a set of complex numbers. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 

CALL Cadd (N, A(*), B(*), Real, Imag) 

• Input Parameters 



N 
A(*) 

B(*) 



Output Parameters 

Real 
Imag 



Number of complex numbers to be added. 

Vector containing the real components of the complex num- 
bers to be added; subscripted from 1 to N. 

Vector containing the imaginary components of the com- 
plex numbers to be added; subscripted from 1 to N. 



Real component of the sum of the complex numbers. 
Imaginary component of the sum of the complex numbers. 



Special Considerations and Programming Hints 

• If there is a wide range in the magnitude of the numbers to be added, smaller quantities 
should be stored first in the arrays A(*) and B(*) to avoid roundoff error. 
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Multiplication with Complex Numbers 

Description 

Subprogram Cmult multiplies two complex numbers. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 

CALL Cmult (Al, Bl, A2, B2, R, I) 

• Input Parameters 



Al, A2 
Bl, B2 

Output Parameters 

R 
1 



Real components of the two complex numbers to be multi- 
plied. 

Imaginary components of the two complex numbers to be 
multiplied. 



Real component of the product. 
Imaginary component of the product. 
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Division with Complex Numbers 

Description 

Subprogram Cdivid divides one complex number by another. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 

CALL Cdivid (Al, Bl, A2, B2, R, I) 

• Input Parameters 



Al 
Bl 
A2 
B2 

Output Parameters 

R 



Real component of the dividend. 
Imaginary component of the dividend. 
Real component of the divisor. 
Imaginary component of the divisor. 



Real component of the quotient. 
Imaginary component of the quotient. 
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Square Root of a Complex Number 

Description 

Subprogram Csqrt finds the square root of a complex number. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 
CALL Csqrt(A, B, R, I) 

• Input Parameters 



A 
B 

• Output Parameters 

R 

I 

Method and Formula 



Real component of the complex number. 
Imaginary component of the complex number. 



Real component of the square root. 
Imaginary component of the square root. 



_ A 2 + B 2 + A A 2 + B 2 - A 

K = ~ ; 1 — ^ 



Specicil Considerations 

• To keep the radical real, the positive sign will always be used with the |a + ib|. 
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Exponential Value of a Complex Number 

Description 

Subprogram Cexp finds the exponential value of a complex number. 

Program Usage 
Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 
CALL Cexp(A, B, R, I) 

• Input Parameters 



A 
B 

• Output Parameters 

R 
1 



Real component of the complex number. 
Imaginary component of the complex number. 



Real component of the exponential. 
Imaginary component of the exponential. 



Method and Formula 

R = e A *cos(B); I = e A *sin(B) 

Special Considerations 

• All of the trigonometric functions require that the argument of the function be given in 
radians. 
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Natural Logarithm of a Complex Number 

Description 

Subprogram Clog finds the natural logarithm of a complex number. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 
CALL Clog (A, B, R, I) 

• Input Parameters 



A 
B 

• Output Parameters 

R 

I 

Method and Formula 

R = LOG( VA 2 + B 2 ); I 



Real component of the complex number. 
Imaginary component of the complex number. 



Real component of the natural logarithm. 
Imaginary component of the natural logarithm. 



= atn (B/A) if A > 

= atn (B/A) + acs ( - 

= atn (B/A) - acs ( - 

= asn(SGN(B))ifA 



1) if A <0andB^0 
1) if A<0andB<0 
= 



Special Considerations 

• All the trigonometric functions require that the argument of the function be given in 
radians. 

• The logarithm function is a multivalued function. When this function is used, the princi- 
pal value of the function is the calculated result. Thus, the function LOG(EXP(Z)) will not 
necessarily return the value Z. Returned phase I is defined on the interval [ - it, + it]. 
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Absolute Value of a Complex Number 

Description 

Subprogram Cabs finds the absolute value of a complex number. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 

CALL Cabs (X, Y, Cabs) 

• Input Parameters 

X Real component of the complex number. 

Y Imaginary component of the complex number. 

• Output Parameter 

Cabs Absolute value of the complex number x + iy. 
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Inverse of a Complex Number 

Description 

Subprogram Cinv finds the inverse of a complex number. 

Program Usage 

Subprogram Utilization 

• File Name - contained in "Complx", disc 2 

• Calling Syntax 
CALL Cinv (X, Y, R, I) 

• Input Parameters 



X 
Y 

• Output Parameters 

R 



Real component of the complex number. 
Imaginary component of the complex number. 



Real component of the inverse. 
Imaginary component of the inverse. 



Method and Formula 

R == X/(X 2 + Y 2 ); I - - Y/(X 2 + Y 2 ) 
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Cosine of a Complex Number 

Description 

Subprogram Ccos finds the cosine of a complex number. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 
CALL Ccos (X, Y, R, I) 

• Input Parameters 

X Real component of the complex number 

Y Imaginary component of the complex number 

• Output Parameters 

R Real component of the cosine. 

I Imaginary component of the cosine. 

Method and Formula 

R = cos (X) cosh (Y); I = - sin (X) sinh (Y) 

Special Considerations 

• All of the trigonometric functions require that the argument of the function be given in 
radians. 

• The cosecant function (R + i * I) = esc (X + i * Y) can be computed by using 
CALL Ccos (X, Y, R, I) 

CALL Cinv (R, I, R, I) 
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Sine of a Complex Number 

Description 

Subprogram Csin finds the sine of a complex number. 
Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 
CALL Csin (X, Y, R, I) 

• Input Parameters 

X Real component of the complex number. 

Y Imaginary component of the complex number. 

• Output Parameters 

R Real component of the sine. 

I Imaginary component of the sine. 

Method and Formula 

R = sin (X) cosh (Y); I = cos (X) sinh (Y) 

Special Considerations 

• All of the trigonometric functions require that the argument of the function be given in 
radians. 

• The secant function (R + i * I) = sec(X + i * Y) can be computed by using 
CALL Csin (X, Y, R, I) 

CALL Cinv (R, I, R, I) 
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Tangent of a Complex Number 

Description 

Subprogram Ctan finds the tangent of a complex number. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 
CALL Ctan (X, Y, R, I) 

• Input parameters 

X Real component of the complex number 

Y Imaginary component of the complex number. 

• Output Parameters 

R Real component of the tangent. 

I Imaginary component of the tangent 

Method and Formula 

Z = (R + i * I) = sin (X + i * Y)/cos (X + i * Y) 

Special Considerations 

• All of the trigonometric functions require that the argument of the function be given in 
radians. 

• The cotangent function (R + i * I) = cot (X + i * Y) can be computed by using 
CALL Ctan (X, Y, R, I) 

CALL Cinv (R, I, R, I) 
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Hyperbolic Cosine of a Complex Number 

Description 

Subprogram Ccosh finds the hyperbolic cosine of a complex number. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 

CALL Ccosh (X, Y, R, I) 

• Input Parameters 

X Real component of the complex number 

Y Imaginary component of the complex number 

• Output Parameters 

R Real component of the hyperbolic cosine 

I Imaginary component of the hyperbolic cosine 

Method and Formula 

R = cosh (X) * cos (Y); I = sinh (X) * sin (Y) 

Special Consideration 

• All of the trigonometric functions require that the argument of the function be given in 
radians. 

• The hyperbolic cosecant (R + i * I) = csch (X + i * Y) can be computed by using 
CALL Ccosh (X, Y, R, I) 

CALL Cinv (R, I, R, I) 
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Hyperbolic Sine of a Complex Number 

Description 

Subprogram Csinh finds the hyperbolic sine of a complex number. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 

CALL Csinh (X, Y, R, 1) 

• Input Parameters 

X Real component of the complex number. 

Y Imaginary component of the complex number. 

• Output Parameters 

R Real component of the hyperbolic sine 

I Imaginary component of the hyperbolic sine 

Method and Formula 

R - sinh (X) * cos (Y); R = cosh (X) * sin (Y) 

Special Consideration 

• All of the trigonometric functions require that the argument of the function be given in 
radians. 

• The hyperbolic secant (R + i * Y) = sech (X + i * Y) can be computed by using 
CALL Csinh (X, Y, R, I) 

CALL Cinv (R, I, R, I) 
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Hyperbolic Tangent of a Complex Number 

Description 

Subprogram Ctanh finds the hyperbolic tangent of a complex number. 

Program Usage 

Subprogram Utilization 

• File Name - contained in "Complx", disc 2 

• Calling Syntax 

CALL Ctanh (X, Y, R, I) 

• Input Parameters 

X Real component of the complex number. 

Y Imaginary component of the complex number. 

• Output Parameters 

R Real component of the hyperbolic tangent. 

I Imaginary component of the hyperbolic tangent. 

Method and Formula 

Z = (R + i * I) = sinh (X + i * Y)/cosh (X + i * Y) 

Special Considerations 

• All of the trigonometric functions require that the argument of the function be given in 
radians. 

• The hyperbolic cotangent (R + i * I) = coth (X + i * Y) can be computed by using 
CALL Ctanh (X, Y, R, I) 

CALL Cinv (R, I, R, I) 
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Rectangular to Polar Conversion 

Description 

Subprogram Rec_pol finds the magnitude and the phase (in the angular unit of the calling 
context) of a complex number or point in rectangular coordinates. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 

CALL (A, B, Magn, Phase) 

• Input Parameters 

A Real component of the complex number, or X-coordinate of 

the point. 

B Imaginary component of the complex number, or Y- 

coordinate of the point. 

• Output Parameters 

Magn Magnitude (orAbsolute Value) of the complex number, or 

p-coordinate of the point. 

Phase Phase of the complex number, or O-coordinate of the point, 

in angular unit of the calling context. 

Method and Formu la 

Magn = VA 2 + B 2 ; Phase* = atn (B/A) if A > 

= atn (B/A) + acs(-l) if A < 
= asn(SGN(B)) if A = 
(* module 2 * tt definition) 

Special Considerations 

• Unlike other complex functions, the angular unit mode can be chosen. 

• Returned phase is defined on the interval [ + 180°, + 180°] or [ - tt, + it] 
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Polar to Rectangular Conversion 

Description 

Subprogram PoLrec finds the complex number (or the rectangular coordinates) of a point 
defined by its magnitude and phase (or polar coordinates) in the angular unit of the calling 
context. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 

CALL PoLrec (A, B, Magn, Phase) 

• Input Parameters 

Magn Magnitude of the complex number, or p-coordinate of the 

point. 

Phase Phase of the complex number, or 6-coordinate of the point 

in angular unit of the calling context. 

• Output Parameters 

A Real component of the complex number, or X coordinate of 

the point. 

B Imaginary component of the complex number or Y coordin- 

ate of the point. 

Method and Formula 

A = Magn * cos (Phase); B = Magn * sin (Phase) 

Special Considerations 

• Unlike other complex functions, the angular unit mode can be chosen. 
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Evaluation of a Complex Polynomial 

Description 

Subprogram Polyev evaluates the complex polynomial 

(a + b i) + (a x + b^Z + ... + (a n + b n i)Z n at a complex point. 

Program Usage 

Subprogram Utilization 

• File Name - contained in file "Complx", disc 2 

• Calling Syntax 

CALL Polyev (N, A, B, Rcoef(*), Icoef(*), Rval, Ival) 

• Input Parameters 



N 
A 

B 

Rcoef(*) 

Icoef(*) 



Degree of the complex polynomial. 

Real component of the complex number at which the 
polynomial is to be evaluated. 

Imaginary component of the complex number at which the 
polynomial is to be evaluated. 

Vector containing the real components of the coefficients of 
the polynomial (Rcoef(O) +Icoef(0)*I) + (Rcoef(l) 
+ Icoef(l)*I)*Z + ... + (Rcoef(N) + Icoef(N)*I)*Z | N; 
subscripted from to N. 

Vector containing the real components of the coefficients of 
the polynomial (Rcoef(O) + Icoef(0)*I) + (Rcoef(l) 
+ Icoef(l)*I)*Z + ... + (Rcoef(N) + Icoef(N)*I)*Z | N; 
subscripted from to N. 



• Output Parameters 

Rval 
Ival 



Real component of the evaluated polynomial. 
Imaginary component of the evaluated polynomial. 
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Fourier Analysis 



Introduction 

Description 

This section contains routines for computing Fourier series coefficients using equally or une- 
qually spaced data, as well as a fast Fourier transform routine. 

Routines 

• Foureq - calculates the Fourier series coefficients using equally-spaced data points. 

• Fourun - calculates the Fourier series coefficients using unequally-spaced data points. 

• Fft - calculates a Fast Fourier Transform from a set of time domain points to a set of 
frequency domain points or vice versa. 
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Fourier Series Coefficients 
for Equally-spaced Data Points 

Description 

This subprogram calculates the Fourier Series Coefficients a, and bj of the Fourier Series 
corresponding to a function f(x) which is specified by n discrete equally-spaced data points 
(x,, y,), i = 1, ..., n. 

Program Usage 

Driver Utilization 

• File Name - "FOUREQ", disc 2 

The driver "FOUREQ" sets up the necessary input parameters for the subprogram Foureq 
and prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Foureq", disc 2 

• Calling Syntax 

CALL Foureq (N, Npts, Init, Incre, Y(*), A(*), B(*)) 

• Input Parameters 

N Highest numbered Fourier Series coefficient. 

Npts Number of data points (must be odd). 

Init Initial domain value Xj. 

Incre Increment between domain values; Ax = x i + 1 — x r 

Y(*) Vector containing range values for data points; subscripted 

from 1 to Npts. 

• Output Parameters 

A(*) Vector containing the A R Fourier Series coefficients; sub- 

scripted from to N. 

B(*) Vector containing the B R Fourier Series coefficients; sub- 

scripted from 1 to N. 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Foureq." 
N = Npts = Incre = 

The data may be corrected from the keyboard (e.g., N = 5, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• Since a Simpson's Method is used in the computation of the coefficients, an odd number 
of data points must be provided. If this is not the case, the following error message is 
printed and the program pauses: 

"ERROR IN SUBPROGRAM Foureq." 

"ODD NUMBER OF DATA POINTS REQUIRED." 

Npts = 

The data may be corrected from the keyboard (e.g., Npts = 15, EXECUTE). When 
CONTINUE is pressed, the program will resume execution at the next line. 

• For valid results, you should provide at least n data points to compute n coefficients. 

• Y(*), A(*) and B(*) must all be dimensioned in the calling program. 

Methods and Formulae 

If g(x) = f(x) cos ( j X j and h(x) = f(x) sin ( ™ x J then: 

ai = 2A3L. { g ( Xl ) + 4g(x 2 ) + 2g(x 3 ) + 4g(x 4 ) + ... + 4g(x n _!) + g(x n )} 

b, - 2Ajl. {h ( Xl ) + 4h(x 2 ) + 2h(x 3 ) + 4h(x 4 ) + ... + 4h(x n _i) + h(x n )} 

Sine and cosine functional values are computed recursively with the following formula: 
,(^„ + !„ - sln (^)cos(^) + cos^n^j) 



sin 



cos 



(^ (J + „, . cos ( ^ ) cos (?a_ a) _ sln (^) sln (^ J ) 
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User Instructions 

1. Make sure Numerical Analysis flexible disc 2 is inserted correctly into the disc drive. Load 
file "FOUREQ" and press RUN. 

2. You will be asked to supply entries for the following items: 

- number of data points to be used in the calculation (must be an odd number) 

- highest coefficient desired in the Fourier Series (For valid results this should be less 
than or equal to the number of data points. ) 

- initial domain value (first x-value) 

- increment between x-values 

- range values (y-values) for each data point. 
Press CONTINUE after each entry. 

3. The program will calculate and print the Fourier Series Coefficients using the given data 
points. 



Fourier Analysis 8-5 



Example 

User entries: 



# OF DATA POINTS^ 3i 
HIGHEST FOURIER SERIES COEI 
INITIAL DOMAIN YALUE- 
INCREMENT^ .1 
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Results: 
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Fourier Series Coefficients 
for Unequally-spaced Data Points 

Description 

This subprogram calculates the Fourier Series Coefficients for a function defined by discrete 
data points (x b y t ), i = l,...,n. The data pairs must be entered such that the x s are discrete, but 
not necessarily equally-spaced, and Xj < Xj + : for i = l,...,n — 1. 

Program Usage 

Driver Utilization 

• File Name - "FOURUN", disc 2 

The driver "FOURUN" sets up the necessary input parameters for the subprogram Fourun 
and prints out the resulting outputs. 

Subprogram Utilization 

• File Name - "Fourun", disc 2 

• Calling Syntax 

CALL Fourun (N, Npts, X(*), Y(*), A(*), B(*)) 

• Input Parameters 

N Highest numbered Fourier Series coefficient. 

Npts Number of data points. 

X(*) Vector containing domain values for data points; sub- 

scripted from 1 to Npts. 

Y(*) Vector containing range values for data points; subscripted 

from 1 to Npts. 

• Output Parameters 

A(*) Vector containing the A k Fourier Series coefficients; sub- 

scripted from to N. 

B(*) Vector containing the B k Fourier Series coefficients; sub- 

scripted from 1 to N. 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Fourun." 

N = Npts - 

X(l) = X(Npts) = 

The data may be corrected from the keyboard (e.g., N = 5, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 

• For valid results, you should provide at least n data points to compute n coefficients. 

• X(*), Y(*), A(*), B(*) must all be dimensioned in the calling subprogram. 

• The data points (x,, y,) should be entered so that the x 5 are discrete and x ( < x, + x for 
i = 1 to n - 1. The points need not be equally spaced. 

Methods and Formulae 

This subprogram calculates the Fourier Series coefficients a s and b, of the Fourier Series cor- 
responding to a function f(x) which is specified by n discrete data points (x,, y { ), i = 1 to n. 

The finite Fourier Series is given by the formula: 

n 

-^2- + S (a;Cos -^- + biSin-i^- ) where the Fourier coefficients a t and b; are: 
2 i-i T T 

a- t = ~y J"* 1 f(x)cos — j — dx for i =0 ton, and 

b| = -j- f*\ f(x)sin j dx for i = 1 to n. 

T specifies the period equivalent to (Xj - x x ) and n indicates the number of coefficients de- 
sired. The coefficients are evaluated by numerically integrating a parabola passing through 
three successive points. Execution time depends on the number of coefficients calculated. 

Formulae: 

k-l 

A, = £ Si for j = 1 to n 



k-l 

B; = 2 T, for j = 1 to n 



8-8 Fourier Analysis 

where: 

(x 1 + 1 -x i42 ) y>-(xi-x i + 2 ) y i+ i + (Xj-x i + 1 ) y i + 2 



Q, 



(X, f i-X i + 2 ) (Xi~X i + 2 ) (Xj-Xi_ 



A = \ (Q, + Q,_j) 



R = -(x iM 2 -x i + 2 2 ) ^-(x^-Xj+g 2 ) yj + t + jx^-x^! 2 ) y i + 2 
(x i + 1 -x i + 2 ) (Xi-X i + 2 ) (Xi-X i + 1 ) 



B = \ (R, + R^) 



(2Ax i , 1 + B)cos[2ir (^pzV)] - (2A Xj + B) cos [2tt (tt^^-) ] 

C _ \ x k x l / \ X k — X t / 

V x k - Xj / 

[2tt l^xTJ ] " 



/ 2ttJ \ 

\X k - Xj/ 

( — V 

\X k - X!/ 



sin 



/ 2itJ y 

\X k - Xj/ 



2A 

— sin 



in ^(^)] 



T, = 



(2Ax i + 1 + B)sin[2Tifrc (~^) ]-(2Ax, + b) sin[2Trfrc (~^) ] 

k (_j_v k ~^~ 

\ X k - Xi / 
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; v cos [ 2lT I v _ v j] + 

( 2ttJ y VXk Xi/ 

\ X W - Xi / 



\ X k - Xj / 



where J is the number of coefficient being computed and frc () indicates fractional part. 



k = 



A = Hi 



1=2 



tt A( x, + 1 3 - x, 3 ) B(x i + 1 2 - Xi 2 ) 

Uj = o I" 7p + MX i + i — X i( 

where A, B are defined above and: 



C = \ (P, + Pi 



a 



(x i + 1 - x i + 2 )x i + 1 x i + 2 Vj - (Xj - Xi + 2 )x i Xi + 2 y 1 + 1 + (x t - x i + 1 )Xj x i + 1 y i + 2 

1 ~ (X i + 1 - X i + 2 ) (^ - X i + 2 ) (Xj - Xj + : ) 

Sine and cosine functional values are computed recursively with the following formulae: 
/ 2ttXj \ / 2ttk, \ ( 2^ \ / 2mc, \ . / 2ttx, \ 



cos 



where i is the data point number and J is the coefficient number. 
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References 

1. Hewlett-Packard 9820A Math Pac, pp. 43-50. 

2. Hamming, R.W., Numerical Methods for Scientists and Engineers, McGraw-Hill, 1962, 
pp. 67-80. 

3. Acton, Forman S., Numerical Methods that Work, Harper and Row, 1970, pp. 221- 
257. 

User Instructions 

1. Make sure Numerical Analysis flexible disc 2 is inserted correctly into the disc drive. Load 
file "FOURUN" and press RUN. 

2. You will be asked to supply entries for the following items: 

- number of data points to be used in the calculation 

- highest Fourier Series Coefficient desired (For valid results, this should be less 
than or equal to the number of data points provided. ) 

- discrete x-values (in increasing order) and y-values for each data point. 
Press CONTINUE after each entry. 

3. The program will print the resulting A k and B k Fourier Series Coefficients. 
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Example 

User entries: 



* OF DATA POINTS^ 37 

HIGHEST FOURIER SERIES COEFFICIENT^ iO 



DATA 


VAU. 


EB 










X ( 


i ) = 





OOOOOOE+00 


Y ( i ) ™ 


i 


50 00 OE+0 


X ( 


2)== 


i 


000 OE-Oi 


Y( 2) = 


c> 


OOOOOOE+0 


X ( 


?\ ') : " 


i 


50 00 OE-Oi 


Y ( 3 ) - 


3 


OOOOOOE+0 


X ( 


4 ) ::: 




00000 OE-Oi 


Y ( 4 ) - 


•■/. 


OOOOOOE+0 


X ( 


5)" : 


''j 
t... 


50 00 0E-0i 


Y ( 3 ) = 


3 


OOOOOOE+0 


X ( 


6 ) - 


3 


00 00 OE-Oi 


Y ( 6 ) = 


"?. 


OOOOOOE+0 


X i - 


7 ) ::: 


4 


00 00 0E-0i 


Y ( 7 ) :•- 


3 


OOOOOOE+0 


X( 


8 ) ::: 


5 


000 OE-Oi 


Y ( 8 ) = 


3 


OOOOOOE+00 


X( 


9):r 


6 


00 OE-Oi 


Y ( 9 ) - 


3 


OOOOOOE+0 


X ( 


1 ) '= 


'? 


00 0E-0i 


Y ( i ) = 


3 


OOOOOOE+0 


X( 


i i ) -- 


8 


00000 OE-Oi 


Y (', i ).■= 


3 


OOOOOOE+00 


';< ( 


1 2 ) - 


9 


000 OE-01 


Y ( i 2 ) - 


3 


OOOOOOE+0 


x i 


't /', ) ::: 


i 


OOOOOOE+00 


Y ( i 3 ) = 


i 


50 00 OE+0 


X( 


1 4 ) ::: 


i 


iO000OE+O0 


Y ( i 4 ) ™ 





OOOOOOE+00 


X ( 


ir>)« 


• i 


2 00 OE + 


Y(i.S) = 





E + 


x>: 


16)"- 


i 


30 00 OE+0 


Y ( i 6 ) ■- 





OOOOOOE+0 


X< 


i 7 ) :: 


:i. 


SO 00 OE+0 


Y ( i 7 ) ~ 





OOOOOOE+0 


X ( 


1 8 ) : '- 


i 


7 E + 


Y ( i 8 > : - 





OOOOOOE+0 


X ( 


1 9 ) - 


i 


90 00 OE+0 


Y ( i 9 ) == 





00 OOF*- 


X ( 


20>~ 


••:> 
/... 


OOOOOOE+00 


Y ( 2 ) - 





OOOOOOE+00 


X ( 


2 1 ) = 


id 


50 00 OE+0 


Y ( 2 i ) :::: 





OOOOOOE+0 


X 1 : 


22 ) ::: 


o 


730 00 OE+0 


Y<22)-- 





i"j o n 1: + 


X ( 


23 X : 


-.^ 


00 0E + 


Y (23)- 





00 OE + 00 


V ( 


24 ) :: 


: 3 


SO 00 OF *fl!i 


Y( 24)-" 





n i) o o o i : mj o 


X i 


■'3 )"• 


",V 


7 E + 


Y \ ", r ■ ) - 


!) 


1 1;:. + 


A ( 


."■;,;_, ) -: 


- 3 


9 ') F * u 


Y ( 2 6 ) - 





E + 


x c 


2/ > 


: 4 


E + 


Y<27) :; ^ 





OOOOOOE+0 


y ( 


•. : 'r ) 

2 ':■'' ) ' r 


'"I 

■ 4 


, ''-OOO 01 Mi ;.i 
"•'■', u (i ! M) 
3 OOu E + I 


Y>;28)- 

Y ( 29 ) ■= 

Y c 3 ) ™ 


[1 





i:" + 
! : + 
OE+00 


X ( 


3? ) ::: 




E + 
2 (j E + 


Y ( 3 i ) -- 
Y ( 32 ) - 





E + 

no oo o Qi.. + o o 


X ( 


33 ) -■■■ 


c 


4 E + 


Y ( 33 ) - 





E + 


X ( 


'■\ .:'). j :: 


c; 


5 :; S3) | M'i '■ I 


Y(34>-~ 





!::'. + 


X ( 
x ( 
X \ 


J'\ I 1 '-, ) " 


- 5 
v." 


830 00 OE+00 
9 E •* 
OOU 01 - (i 


Y(3S>= 

Y ( 3 6 ) :::: 

Y ( 37 ) ■■-■ 




1) 


OOOOOOE+0 
OE + (: 
50 00 OE+0 
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Results: 



A COEEF ICIENT 



!< COEFFICIENT 



'"■ I) 'J l.i; "U 1 

:■:: 26 2 33 BE ■•■■01. 
4 . i. 1 97(:vE-0i 
■\ ?>fi- / -I '■>" ■} (■. 

2 n: ; ;;vi.3SE--u.i. 

i . 6i4 37F ■•■ 01 
2 46 ; ;-94E- :i.2 
:'i. i 7 ' : ° 1 E ii :i. 
v.63:i. 3377: ■ 2 
2 7773250E-1.7 
'■''. 36 777 IE- 02 



•<{ / ' ; i ■ 


i o (.; , 


- i 

- (! :i 


b .307: 


3S0F 


- :). 


3 . '".- MV 


.5 4 E 


- :i 


.} . -v .,. .■■ { 

i . 6 6 6 ■■ 

:i. . ° .'; ■■:: i 
1 . 2' 75/ 


'. •:> 1 ■;:. 

'-23 
i '::> '".-' r:' 

>92E 


■" U 7 
•314 
... ()fl 

- :t 
•01 
■ 1 
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Fast Fourier Transform 

Description 

This subprogram calculates a Fast Fourier Transform from a set of time domain points to a set 
of frequency domain points or the inverse Fast Fourier Transform, calculating the set of time 
domain points from a set of frequency domain points. 

Program Usage 

Driver Utilization 

• File Name - "FFT", disc 2 

The driver "FFT" sets up the necessary input parameters for the subprogram Fft and prints 
out the resulting outputs. 

Subprogram Utilization 

• File Name -"Fft", disc 2 

• Calling Syntax 

CALL Fft(N, Power, Fig, R(*), I(*)) 

• Input Parameters 



N 

Power 
Fig 



R(*Y 



K*) 1 



Output Parameters 
RW 1 
K*) 1 



(number of time domain points)/2; or (number of frequency 
coefficient pairs)/2 + 1. 
M _ 2( p ° wer ~i> 

If Fig = 1, calculate Fast Fourier Transform from a set of 
time domain points to a set of frequency domain points. 
If Fig = - 1, calculate Inverse Fast Fourier Transform from 
a set of frequency domain points to a set of time domain 
points. 

Vector subscripted from 1 to N; as time domain data, R(*) 
contains the odd-numbered data points, i.e., Rj = time do- 
main in data for i = 1 to n; as frequency domain data, Rj 
contains the DC term, and Rj = real component of data 
(j - l)forj = 2 ton. 

Vector subscripted from 1 to N; as time domain data, I(*) 
contains the even-numbered data points, i.e., Ij = time do- 
main data (2*j) for j = 1 to n; as frequency domain data, I x 
contains the maximum frequency term, and I s = imaginary 
component of data (j - 1) for j = 2 to n. 



Vector subscripted from 1 to N (see R(*) as an input param- 
eter). 

Vector subscripted from 1 to N (see I(*) as an input param- 
eter). 



1 R{*) and I(*) are both input and output parameters. 
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Special Considerations and Programming Hints 

• Upon entry into the subprogram, there is a bad data check. If the subprogram detects 
"nonsense" data, the following error message is printed and the program pauses: 

"ERROR IN SUBPROGRAM Fft." 

N = Fig = Power = 

The data may be corrected from the keyboard (e.g., Fig = 1, EXECUTE). When CON- 
TINUE is pressed, the program will resume execution at the next line. 



Data Points 






(Full Precision Arrays) 


FFT 


IFT 


256 


4.3 sec 


4.0 sec 


512 


9.1 sec 


8.7 sec 


1024 


19.6 sec 


18.6 sec 
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Methods and Formulae 

This method uses a modification of the basic FFT algorithm. The modified algorithm takes 
advantage of the fact that series data will be real and the space normally reserved for the 
imaginary part of the complex sequence can be used to calculate a double-length real trans- 
form. This is represented for two N length transforms as: 

Z(n) = X(n) + i Y(n) =s n < N data points 

The transform is: 

Z(m) = X(m) + iY(m) 
where 

Y , , Z(m) + Z*(N - m) 
X(m) = 2 

Y(m) = -^ p — Z* is the complex conjugate of Z. 

The time series F(n) is given by: 
F(n) - X(2n) + Y (2n + 1) 



The advantage gained from this adaptation of the general FFT algorithm for time series data 
are: 

1. A transform of twice the length can be handled with no increase in storage for input 
data. 

2. Since the calculation of the transform is treated as an interactive process, intermediate 
and final results are stored in the same locations used for input. 

The transform of this is: 



N-l N-l 

X X(2n)w mn + 2 

n=0 n=0 



F(m) = 2 X(2n)w mn + 2 Y(2n + l)w mn 



N-l N-l 

= X X(p)w 2mp + 2 Y(p)w 2mp w n 

p=0 p-0 
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and: 

F(m) = X(m) + w m Y(m) (1) 

F(N - m) = X*(m) - [w m Y(m)]* (2) 
Similarly the inverse transform may be obtained from (1) and (2): 
7 , , F(m) + F*(N - m) . _ m F(m) - F*(N - m) 

Z(N-m)^[ F(m) + F 2 - (N - m) ]* _ iw - m [M^(N_jnl]* 

This is simply an interchange of Z(m) and F(m) in (1) and (2), and substitution of ( - w~ m ) for 
w m . 

Note: 

1. Since F(0) and F(N) are real only, F(N) can be stored in the imaginary location of F(0), 
i.e.,F(l). 

2. w m = c - 2mm/2N _ 7his is half the minimum value of rotation normally used in an N point 
transform. 

3. * = complex conjugate. 

References 

1. Brigham, E.O., The Fast Fourier Transform, Englewood Cliffs, N.J.: Prentice Hall, Inc., 
1974, Chapter 10. 

2. Cooley, J.W., and Tukey, J.W., "An Algorithm for Machine Calculation of Complex 
Fourier Series," Math. Computation, Vol. 19 (April 1965), pp. 297-301. 
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User Instructions 

1. Make sure Numerical Analysis flexible disc 2 is inserted correctly into the disc drive. Load 
file "FFT" and press RUN. 

2. You can append subprogram Fft by answering "no" ("N") when asked, "HAS Fft 
ALREADY BEEN APPENDED (Y/N)?" This will cause file "Fft" to be linked to the end 
of the driver program. This only needs to be done when the driver program (FFT) is 
initially loaded or reloaded. 

3. You will be asked to supply entries for the following items: 
If time data is to be entered: 

- number of time domain points (power of 2 from 4 to 1024) to be used 
in the calculations 

- time domain data points. 

The program will print the DC term, the maximum frequency, and the real and imagin- 
ary components of the resulting frequency data, 
or 
If frequency data is to be entered: 

- value of the number of coefficient pairs times 2 plus 2 (power of 2 
between 4 and 1024) 

- DC term 

- maximum frequency term 

- real and imaginary coefficients for the frequency domain data points. 

The program will print the time domain data points obtained from the inverse Fast 
Fourier Transform. 
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Examples 

• Time to frequency domain 

• Time domain data was generated using the following formula: 



POINT(I) = SIN((I-l)*ir/8), forI = ltol6. 



User entries: 



i i.+ 1'HlA P0INT0- 16 



lAdi POHAiN DATA 



PO i' NT ' l "■> :•■ U U i.l «.» i;-. -*- it 



POINT ( 
POINT i 
POINT( 
POINT ( 
POINT* 
POINT-; 
P 1 N T ( 
POINT ( 
POIN 



9 ) "- 3 . 3 268 34 !:O-0 1 

3) ••"■ 7 07 a 68E-01 

7 > ™ i I) U U P + 

6) ::: - V. 238793E-0 1 

7 )■■■,■ 7 . 710681::.- 01 

3)- 3 876834E- 01 

v "j -■ i.) u it (O) UP t : .i U 

-I- U / ' -.J . i.'.'...'..M..'i.'*TI.„ u .1. 



PO.I.NI ( i 1 )™- 7 07 1 OOiil -0:1. 



•' i ::.: •• '■■' .J ■< )•■! 



I" U .1. N ! ( Id) ■■■■ 

POINT ( 13 0-' 1 . HO 1" + 00 

P J.' N T ( 14) ■■■■■ •••■ 9 2 3 8 7 9B P •- .1. 

POINT*: IS)™ -7 . 071068E- 01 

P O I N T ( 16) ™ ■- 3 . 8 2 6 8 3 4 IP ■■■■ 1 



Results: 



DC TERN™ 

HAX FREQUENCY™ 

FREQUENCY DOMAIN: 

REAL. 

1 2.775358E -17 

. OOOOOOE+0 

3 4 . 274748E-17 

4 0. 00O00OE + 0O 
3 "-9. 8 2 58 6 3E -17 

6 . OOOOOOE + 00 

7 2 . 775358E-17 



! MAO I NARY 

-:i. . OQOOEMi 
. IP + 
2 78397P :■- 9 

OOOOOOE+00 

1 538152E-08 
. OOOOOOE + 00 
2. 7 8626 IE -0 8 



• Time to frequency domain 
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User entries: 



# OF DATA POINTS-- 16 



TIME DOMAIN DATA: 



POINT* 


i ) 


POINT ( 


2) 


POINT ( 


3 ) 


POINT* 


■4) 


POINT ( 


S ) 


POINT < 


6 ) 


POINT* 


"7 "i 
Q 'i 


1" U .1. IN ! '-. 

POINT ( 


9) 


POINTC 


i ) 


1'U.I.NT* 


i i ) 


P I NT ( 


12) 


POINT ( 


13) 


POINT* 


i 4 ) 


POINT* 


15) 


POINT* 


16 ) 



E 
3 6600E 
4 1420 E 
3O66O0E 
E 
4120 0E 
E 
412 0E 
E 
3 6600E 
41420 OE 
3 06600E 
1 . ( j E 
5.412000 E 
. E 
5.412000 E 



:i. 
1 
i 
1 
1 
5 




+ 
•*• 
+ 
+ 
+ 
■■■■ 1 
+ 
•■■■ 1 
+ 
+ 
+ 
+ 
+ 
■- 1 
+ 
■- 1 



Results: 



DC TERM- 

MAX FREQUENCY™ 



FREQUENCY DOMAIN. 
PEA!... 



: MAG I NARY 



i. tt 1 i ■•}■ () 
E + 

. 3 E '<■ 
6 i 344 / 6i: •- 'i 

. o i) n o o n e -i o o 

1 . 5 2 2 ''4' U ! > 



1 . () 1 E 




-i .3394531 

. 
■•6 1344761 

. 



1 



.U'd 



>4l. 



+ 

+ 

■- 6 

+ 

■■•■ 6 

:• U 

.... c; 
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• Frequency to time domain 



User entries: 



# F CO E: F F I C I E N T P A I R S * 2 ) + 2 ■■ 



:i. 6 



DC TERM-" 

MAX FREQUENCY TERM= 



EQUENCY DOMAIN DATA 



REA!...< 
REAL'! 
EEA!...( 
REAL( 
REAL. ( 
REALC 
REAL( 



1 > ■ 
2) ■■ 
3 ) ; 

4): 
5): 

6 ) • 

7 t : ' 



0000 i 
000001: 
33950 01 
00 01 
13450 01 



■♦■ 
+ 
-06 
•t- 
-06 



oooooo e:+ oo 

502200 E- 05 



:i:hag( i )■■■■■■■■ ■■■■! . 0000001;:+ 00 

TMAG< 2)=-= . 000000 E+00 

IHAGC 3>^-i . 33950 0E- 06 

IMAG< 4)^ . 000000EHM) 

IMAG< 5) =-■■•■ 6. i345!)0E -06 

IMACX 6)™ 000000 E+ 00 

J: M A G ( 7 ) "- ■••■ 1.502200 e: ■•■• (j 5 



Results: 



i .1. Ml.' DOMAIN : 



DA 
DA 
DA 

D A 
DA 
D A 
DA 
DA 
DA 
Jin 
PA 
DA 
DA 
DA 
DA 
D A 



! H 
!"A 
I'A 
I'A 

T A 



'!" A 



I'A POINT'' 
l'A i'OLN 

t'a p o:i: m< 

I'A PO '.N'T < 
PC) IN' 

point; 

POINT ( 
POINK 
P O I N f C 
I A P f ) I N f ( 
T'A POINT'; 
POINT ( 
POINT ( 



POIN' 



TA POINT' 
f i-'i i ;> 1 N )' ( 



;j ):::: 9 . 9 9 <: > 8 9 8 E - 1 

2 ) ' ::: 1 . 3 6587EHM) 

3) - 1 . 4 14 186 E* 

4) : - i . 3 65G7E-V-0 

5):-: 9 99939911 -01 

6)"" S.4U945E -01 

7)"- 1 . 110923E - 16 

3/ "" •■•'::> 411945E- 01 
9 ■;, :: ~ ... o o 9 9 o 9 o j::- ... q ., 

10 ) -'• 1 . 3 6587EM)0 
ii ) ^"1 . 4 14186 E* 
12) ■•■ •■•■! . 3 6587E>0 
13)^-9 999 898E-01 
:i.4 ) ■■■■■■■■?, 4119451::.'- 01 
15) --1 1 10 2 231) -16 



.1. 6 ) ■■■■ 



41:s94 c 



t ;-.... n ■; 



m 
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